/src/enzo/po_solver.C
C | 4061 lines | 2455 code | 1027 blank | 579 comment | 239 complexity | 6ecdafb4c4c038fa6c6a3d50b6d6f2ca MD5 | raw file
Possible License(s): BSD-3-Clause, LGPL-2.1, GPL-2.0
Large files files are truncated, but you can click here to view the full file
- /* THIS FILE HAS BEEN AUTO-GENERATED. DO NOT EDIT. */
- /* This is C++ code to read HDF5 files for
- reaction rates, cooling rates, and initial
- conditions for the chemical network defined
- by the user. In addition, this contains
- code for calculating temperature from the
- gas energy and computing the RHS and the
- Jacobian of the system of equations which
- will be fed into the solver.
- */
- #include "po_solver.h"
- po_data *po_setup_data(void) {
- po_data *data = (po_data *) malloc(sizeof(po_data));
- /* Temperature-related pieces */
- data->bounds[0] = 1.0;
- data->bounds[1] = 1e+12;
- data->nbins = 1024 - 1;
- data->dbin = (log(data->bounds[1]) - log(data->bounds[0])) / data->nbins;
- data->idbin = 1.0L / data->dbin;
- /* Redshift-related pieces */
- data->z_bounds[0] = 0.0;
- data->z_bounds[1] = 15.0;
- data->n_zbins = 100 - 1;
- data->d_zbin = (log(data->z_bounds[1] + 1.0) - log(data->z_bounds[0] + 1.0)) / data->n_zbins;
- data->id_zbin = 1.0L / data->d_zbin;
-
- po_read_rate_tables(data);
- fprintf(stderr, "Successfully read in rate tables.\n");
- po_read_cooling_tables(data);
- fprintf(stderr, "Successfully read in cooling rate tables.\n");
- return data;
- }
- int po_main(int argc, char** argv)
- {
- po_data *data = po_setup_data();
- /* Initial conditions */
- hid_t file_id = H5Fopen("po_initial_conditions.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
- if (file_id < 0) {fprintf(stderr, "Failed to open "
- "po_initial_conditions.h5 so dying.\n");
- return(1);}
- /* Allocate the correct number of cells */
- hsize_t dims; /* We have flat versus number of species */
- /* Check gas energy to get the number of cells */
- fprintf(stderr, "Getting dimensionality from ge:\n");
- herr_t status = H5LTget_dataset_info(file_id, "/ge", &dims, NULL, NULL);
- if(status == -1) {
- fprintf(stderr, "Error opening initial conditions file.\n");
- return 1;
- }
- fprintf(stderr, " ncells = % 3i\n", (int) dims);
- data->ncells = dims;
- int N = 16;
- double *atol, *rtol;
- atol = (double *) alloca(N * dims * sizeof(double));
- rtol = (double *) alloca(N * dims * sizeof(double));
- double *tics = (double *) alloca(dims * sizeof(double));
- double *ics = (double *) alloca(dims * N * sizeof(double));
- double *input = (double *) alloca(dims * N * sizeof(double));
-
- unsigned int i = 0, j;
-
- fprintf(stderr, "Reading I.C. for /HI\n");
- H5LTread_dataset_double(file_id, "/HI", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 1.00794; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "HI[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /HII\n");
- H5LTread_dataset_double(file_id, "/HII", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 1.00794; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "HII[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /HeI\n");
- H5LTread_dataset_double(file_id, "/HeI", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 4.002602; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "HeI[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /HeII\n");
- H5LTread_dataset_double(file_id, "/HeII", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 4.002602; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "HeII[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /HeIII\n");
- H5LTread_dataset_double(file_id, "/HeIII", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 4.002602; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "HeIII[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /de\n");
- H5LTread_dataset_double(file_id, "/de", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 1.0; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "de[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /ge\n");
- H5LTread_dataset_double(file_id, "/ge", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 1.0; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "ge[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /OI\n");
- H5LTread_dataset_double(file_id, "/OI", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "OI[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /OII\n");
- H5LTread_dataset_double(file_id, "/OII", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "OII[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /OIII\n");
- H5LTread_dataset_double(file_id, "/OIII", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "OIII[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /OIV\n");
- H5LTread_dataset_double(file_id, "/OIV", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "OIV[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /OV\n");
- H5LTread_dataset_double(file_id, "/OV", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "OV[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /OVI\n");
- H5LTread_dataset_double(file_id, "/OVI", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "OVI[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /OVII\n");
- H5LTread_dataset_double(file_id, "/OVII", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "OVII[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /OVIII\n");
- H5LTread_dataset_double(file_id, "/OVIII", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "OVIII[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- fprintf(stderr, "Reading I.C. for /OIX\n");
- H5LTread_dataset_double(file_id, "/OIX", tics);
- for (j = 0; j < dims; j++) {
- ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
- atol[j * N + i] = tics[j] * 1e-09;
- rtol[j * N + i] = 1e-09;
- if(j==0) {
- fprintf(stderr, "OIX[0] = %0.3g, atol => % 0.16g\n",
- tics[j], atol[j]);
- }
- }
- i++;
-
- H5Fclose(file_id);
- double dtf = 3.1557e13;
- double dt = -1.0;
- double z = -1.0;
- for (i = 0; i < dims * N; i++) input[i] = ics[i];
- double ttot;
- ttot = dengo_evolve_po(dtf, dt, z, input, rtol, atol, dims, data);
- /* Write results to HDF5 file */
- file_id = H5Fcreate("po_solution.h5", H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
- hsize_t dimsarr[1];
- dimsarr[0] = dims;
- i = 0;
-
- double HI[dims];
- for (j = 0; j < dims; j++) {
- HI[j] = input[j * N + i] * 1.00794;
- }
- fprintf(stderr, "Writing solution for /HI\n");
- H5LTmake_dataset_double(file_id, "/HI", 1, dimsarr, HI);
- i++;
-
- double HII[dims];
- for (j = 0; j < dims; j++) {
- HII[j] = input[j * N + i] * 1.00794;
- }
- fprintf(stderr, "Writing solution for /HII\n");
- H5LTmake_dataset_double(file_id, "/HII", 1, dimsarr, HII);
- i++;
-
- double HeI[dims];
- for (j = 0; j < dims; j++) {
- HeI[j] = input[j * N + i] * 4.002602;
- }
- fprintf(stderr, "Writing solution for /HeI\n");
- H5LTmake_dataset_double(file_id, "/HeI", 1, dimsarr, HeI);
- i++;
-
- double HeII[dims];
- for (j = 0; j < dims; j++) {
- HeII[j] = input[j * N + i] * 4.002602;
- }
- fprintf(stderr, "Writing solution for /HeII\n");
- H5LTmake_dataset_double(file_id, "/HeII", 1, dimsarr, HeII);
- i++;
-
- double HeIII[dims];
- for (j = 0; j < dims; j++) {
- HeIII[j] = input[j * N + i] * 4.002602;
- }
- fprintf(stderr, "Writing solution for /HeIII\n");
- H5LTmake_dataset_double(file_id, "/HeIII", 1, dimsarr, HeIII);
- i++;
-
- double de[dims];
- for (j = 0; j < dims; j++) {
- de[j] = input[j * N + i] * 1.0;
- }
- fprintf(stderr, "Writing solution for /de\n");
- H5LTmake_dataset_double(file_id, "/de", 1, dimsarr, de);
- i++;
-
- double ge[dims];
- for (j = 0; j < dims; j++) {
- ge[j] = input[j * N + i] * 1.0;
- }
- fprintf(stderr, "Writing solution for /ge\n");
- H5LTmake_dataset_double(file_id, "/ge", 1, dimsarr, ge);
- i++;
-
- double OI[dims];
- for (j = 0; j < dims; j++) {
- OI[j] = input[j * N + i] * 15.9994;
- }
- fprintf(stderr, "Writing solution for /OI\n");
- H5LTmake_dataset_double(file_id, "/OI", 1, dimsarr, OI);
- i++;
-
- double OII[dims];
- for (j = 0; j < dims; j++) {
- OII[j] = input[j * N + i] * 15.9994;
- }
- fprintf(stderr, "Writing solution for /OII\n");
- H5LTmake_dataset_double(file_id, "/OII", 1, dimsarr, OII);
- i++;
-
- double OIII[dims];
- for (j = 0; j < dims; j++) {
- OIII[j] = input[j * N + i] * 15.9994;
- }
- fprintf(stderr, "Writing solution for /OIII\n");
- H5LTmake_dataset_double(file_id, "/OIII", 1, dimsarr, OIII);
- i++;
-
- double OIV[dims];
- for (j = 0; j < dims; j++) {
- OIV[j] = input[j * N + i] * 15.9994;
- }
- fprintf(stderr, "Writing solution for /OIV\n");
- H5LTmake_dataset_double(file_id, "/OIV", 1, dimsarr, OIV);
- i++;
-
- double OV[dims];
- for (j = 0; j < dims; j++) {
- OV[j] = input[j * N + i] * 15.9994;
- }
- fprintf(stderr, "Writing solution for /OV\n");
- H5LTmake_dataset_double(file_id, "/OV", 1, dimsarr, OV);
- i++;
-
- double OVI[dims];
- for (j = 0; j < dims; j++) {
- OVI[j] = input[j * N + i] * 15.9994;
- }
- fprintf(stderr, "Writing solution for /OVI\n");
- H5LTmake_dataset_double(file_id, "/OVI", 1, dimsarr, OVI);
- i++;
-
- double OVII[dims];
- for (j = 0; j < dims; j++) {
- OVII[j] = input[j * N + i] * 15.9994;
- }
- fprintf(stderr, "Writing solution for /OVII\n");
- H5LTmake_dataset_double(file_id, "/OVII", 1, dimsarr, OVII);
- i++;
-
- double OVIII[dims];
- for (j = 0; j < dims; j++) {
- OVIII[j] = input[j * N + i] * 15.9994;
- }
- fprintf(stderr, "Writing solution for /OVIII\n");
- H5LTmake_dataset_double(file_id, "/OVIII", 1, dimsarr, OVIII);
- i++;
-
- double OIX[dims];
- for (j = 0; j < dims; j++) {
- OIX[j] = input[j * N + i] * 15.9994;
- }
- fprintf(stderr, "Writing solution for /OIX\n");
- H5LTmake_dataset_double(file_id, "/OIX", 1, dimsarr, OIX);
- i++;
-
- double temperature[dims];
- for (j = 0; j < dims; j++) {
- temperature[j] = data->Ts[j];
- }
- H5LTmake_dataset_double(file_id, "/T", 1, dimsarr, temperature);
- double time[1];
- time[0] = ttot;
- double timestep[1];
- timestep[0] = dt;
- H5LTset_attribute_double(file_id, "/", "time", time, 1);
- H5LTset_attribute_double(file_id, "/", "timestep", timestep, 1);
- H5Fclose(file_id);
-
- return 0;
- }
-
- double dengo_evolve_po (double dtf, double &dt, double z, double *input,
- double *rtol, double *atol, long long dims, po_data *data) {
- int i, j;
- hid_t file_id;
- /* fprintf(stderr, " ncells = % 3i\n", (int) dims); */
- int N = 16;
- rhs_f f = calculate_rhs_po;
- jac_f jf = calculate_jacobian_po;
- if (dt < 0) dt = dtf; // / 1e3;
- data->current_z = z;
- int niter = 0;
- int siter = 0;
- double ttot = 0;
- double *scale = (double *) alloca(dims * N * sizeof(double));
- double *prev = (double *) alloca(dims * N * sizeof(double));
- for (i = 0; i < dims * N; i++) scale[i] = input[i];
- for (i = 0; i < dims * N; i++) prev[i] = input[i];
- double *u0 = (double *) alloca(N*dims*sizeof(double));
- double *s = (double *) alloca(N*sizeof(double));
- double *gu = (double *) alloca(N*dims*sizeof(double));
- double *Ju = (double *) alloca(N*N*dims*sizeof(double));
- double floor_value = 1e-25;
- while (ttot < dtf) {
- int rv = BE_chem_solve(f, jf, input, dt, rtol, atol, dims, N,
- scale, (void *) data, u0, s, gu, Ju);
- /*
- fprintf(stderr, "Return value [%d]: %i. %0.5g / %0.5g = %0.5g (%0.5g)\n",
- niter, rv, ttot, dtf, ttot/dtf, dt);
- fprintf(stderr, "Value[0] = %0.5g %0.5g\n",
- input[0], prev[0]);
- */
- for (i = 0; i < dims * N; i++) {
- if (input[i] < 0) {
- rv = 1;
- break;
- }
- }
- if (rv == 0) {
- if (siter == 500000) break;
- siter++;
- if (siter % 100000 == 0) {
- fprintf(stderr, "Successful Iteration[%d]: (%0.4g) %0.16g / %0.16g\n",
- siter, dt, ttot, dtf);
- }
- ttot += dt;
- dt = DMIN(dt * 1.1, dtf - ttot);
-
- for (i = 0; i < dims * N; i++) prev[i] = input[i];
- for (i = 0; i < dims * N; i++) {
- if (input[i] < floor_value) {
- input[i] = floor_value;
- }
- scale[i] = input[i];
- }
- } else {
- dt /= 2.0;
- for (i = 0; i < dims * N; i++) {
- input[i] = prev[i];
- scale[i] = input[i];
- }
- if (dt/dtf < 1e-15) {
- fprintf(stderr, "Dying! dt/dtf = %0.5g\n", dt/dtf);
- break;
- }
- }
- /* if ((dt/dtf < 1e-3) && ((dtf - ttot) != 0.0)) {
- fprintf(stderr, "Relative dt is small! (dt = %0.5g, dtf = %0.5g, dt/dtf = %0.5g)\n", dt, dtf, dt/dtf);
- } */
- niter++;
- }
- /* fprintf(stderr, "End: %0.5g / %0.5g (%0.5g)\n",
- ttot, dtf, dtf-ttot); */
- return ttot;
- }
-
- void po_read_rate_tables(po_data *data)
- {
- hid_t file_id = H5Fopen("po_tables.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
- /* Allocate the correct number of rate tables */
- H5LTread_dataset_double(file_id, "/HeI_i", data->r_HeI_i);
- H5LTread_dataset_double(file_id, "/HeI_pi", data->r_HeI_pi);
- H5LTread_dataset_double(file_id, "/HeII_i", data->r_HeII_i);
- H5LTread_dataset_double(file_id, "/HeII_pi", data->r_HeII_pi);
- H5LTread_dataset_double(file_id, "/HeII_r", data->r_HeII_r);
- H5LTread_dataset_double(file_id, "/HeIII_r", data->r_HeIII_r);
- H5LTread_dataset_double(file_id, "/HI_i", data->r_HI_i);
- H5LTread_dataset_double(file_id, "/HI_pi", data->r_HI_pi);
- H5LTread_dataset_double(file_id, "/HII_r", data->r_HII_r);
- H5LTread_dataset_double(file_id, "/OI_i", data->r_OI_i);
- H5LTread_dataset_double(file_id, "/OI_pi", data->r_OI_pi);
- H5LTread_dataset_double(file_id, "/OII_i", data->r_OII_i);
- H5LTread_dataset_double(file_id, "/OII_pi", data->r_OII_pi);
- H5LTread_dataset_double(file_id, "/OII_r", data->r_OII_r);
- H5LTread_dataset_double(file_id, "/OIII_i", data->r_OIII_i);
- H5LTread_dataset_double(file_id, "/OIII_pi", data->r_OIII_pi);
- H5LTread_dataset_double(file_id, "/OIII_r", data->r_OIII_r);
- H5LTread_dataset_double(file_id, "/OIV_i", data->r_OIV_i);
- H5LTread_dataset_double(file_id, "/OIV_pi", data->r_OIV_pi);
- H5LTread_dataset_double(file_id, "/OIV_r", data->r_OIV_r);
- H5LTread_dataset_double(file_id, "/OIX_r", data->r_OIX_r);
- H5LTread_dataset_double(file_id, "/OV_i", data->r_OV_i);
- H5LTread_dataset_double(file_id, "/OV_pi", data->r_OV_pi);
- H5LTread_dataset_double(file_id, "/OV_r", data->r_OV_r);
- H5LTread_dataset_double(file_id, "/OVI_i", data->r_OVI_i);
- H5LTread_dataset_double(file_id, "/OVI_pi", data->r_OVI_pi);
- H5LTread_dataset_double(file_id, "/OVI_r", data->r_OVI_r);
- H5LTread_dataset_double(file_id, "/OVII_i", data->r_OVII_i);
- H5LTread_dataset_double(file_id, "/OVII_pi", data->r_OVII_pi);
- H5LTread_dataset_double(file_id, "/OVII_r", data->r_OVII_r);
- H5LTread_dataset_double(file_id, "/OVIII_i", data->r_OVIII_i);
- H5LTread_dataset_double(file_id, "/OVIII_pi", data->r_OVIII_pi);
- H5LTread_dataset_double(file_id, "/OVIII_r", data->r_OVIII_r);
- H5Fclose(file_id);
- }
- void po_read_cooling_tables(po_data *data)
- {
- hid_t file_id = H5Fopen("po_tables.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
- /* Allocate the correct number of rate tables */
- H5LTread_dataset_double(file_id, "/compton_comp",
- data->c_compton_comp);
- H5LTread_dataset_double(file_id, "/HeI_c_HeI_c",
- data->c_HeI_c_HeI_c);
- H5LTread_dataset_double(file_id, "/HeI_ph_HeI_ph",
- data->c_HeI_ph_HeI_ph);
- H5LTread_dataset_double(file_id, "/HeII_c_HeII_c",
- data->c_HeII_c_HeII_c);
- H5LTread_dataset_double(file_id, "/HeII_ph_HeII_ph",
- data->c_HeII_ph_HeII_ph);
- H5LTread_dataset_double(file_id, "/HeIII_c_HeIII_c",
- data->c_HeIII_c_HeIII_c);
- H5LTread_dataset_double(file_id, "/HI_c_HI_c",
- data->c_HI_c_HI_c);
- H5LTread_dataset_double(file_id, "/HI_ph_HI_ph",
- data->c_HI_ph_HI_ph);
- H5LTread_dataset_double(file_id, "/HII_c_HII_c",
- data->c_HII_c_HII_c);
- H5LTread_dataset_double(file_id, "/OI_c_OI_c",
- data->c_OI_c_OI_c);
- H5LTread_dataset_double(file_id, "/OI_ph_OI_ph",
- data->c_OI_ph_OI_ph);
- H5LTread_dataset_double(file_id, "/OII_c_OII_c",
- data->c_OII_c_OII_c);
- H5LTread_dataset_double(file_id, "/OII_ph_OII_ph",
- data->c_OII_ph_OII_ph);
- H5LTread_dataset_double(file_id, "/OIII_c_OIII_c",
- data->c_OIII_c_OIII_c);
- H5LTread_dataset_double(file_id, "/OIII_ph_OIII_ph",
- data->c_OIII_ph_OIII_ph);
- H5LTread_dataset_double(file_id, "/OIV_c_OIV_c",
- data->c_OIV_c_OIV_c);
- H5LTread_dataset_double(file_id, "/OIV_ph_OIV_ph",
- data->c_OIV_ph_OIV_ph);
- H5LTread_dataset_double(file_id, "/OIX_c_OIX_c",
- data->c_OIX_c_OIX_c);
- H5LTread_dataset_double(file_id, "/OV_c_OV_c",
- data->c_OV_c_OV_c);
- H5LTread_dataset_double(file_id, "/OV_ph_OV_ph",
- data->c_OV_ph_OV_ph);
- H5LTread_dataset_double(file_id, "/OVI_c_OVI_c",
- data->c_OVI_c_OVI_c);
- H5LTread_dataset_double(file_id, "/OVI_ph_OVI_ph",
- data->c_OVI_ph_OVI_ph);
- H5LTread_dataset_double(file_id, "/OVII_c_OVII_c",
- data->c_OVII_c_OVII_c);
- H5LTread_dataset_double(file_id, "/OVII_ph_OVII_ph",
- data->c_OVII_ph_OVII_ph);
- H5LTread_dataset_double(file_id, "/OVIII_c_OVIII_c",
- data->c_OVIII_c_OVIII_c);
- H5LTread_dataset_double(file_id, "/OVIII_ph_OVIII_ph",
- data->c_OVIII_ph_OVIII_ph);
- H5Fclose(file_id);
- }
-
- void po_calculate_temperature(po_data *data,
- double *input, int nstrip, int nchem)
- {
- int i, j;
- double density;
- double kb = 1.3806504e-16; // Boltzmann constant [erg/K]
- double mh = 1.67e-24;
- double gamma = 5.e0/3.e0;
-
- /* Calculate total density */
- double HI;
- double HII;
- double HeI;
- double HeII;
- double HeIII;
- double de;
- double ge;
- double OI;
- double OII;
- double OIII;
- double OIV;
- double OV;
- double OVI;
- double OVII;
- double OVIII;
- double OIX;
- for (i = 0; i<nstrip; i++) {
- j = i * nchem;
- HI = input[j];
- /*fprintf(stderr, "HI[%d] = % 0.16g\n",
- i, HI);*/
- j++;
-
- HII = input[j];
- /*fprintf(stderr, "HII[%d] = % 0.16g\n",
- i, HII);*/
- j++;
-
- HeI = input[j];
- /*fprintf(stderr, "HeI[%d] = % 0.16g\n",
- i, HeI);*/
- j++;
-
- HeII = input[j];
- /*fprintf(stderr, "HeII[%d] = % 0.16g\n",
- i, HeII);*/
- j++;
-
- HeIII = input[j];
- /*fprintf(stderr, "HeIII[%d] = % 0.16g\n",
- i, HeIII);*/
- j++;
-
- de = input[j];
- /*fprintf(stderr, "de[%d] = % 0.16g\n",
- i, de);*/
- j++;
-
- ge = input[j];
- /*fprintf(stderr, "ge[%d] = % 0.16g\n",
- i, ge);*/
- j++;
-
- OI = input[j];
- /*fprintf(stderr, "OI[%d] = % 0.16g\n",
- i, OI);*/
- j++;
-
- OII = input[j];
- /*fprintf(stderr, "OII[%d] = % 0.16g\n",
- i, OII);*/
- j++;
-
- OIII = input[j];
- /*fprintf(stderr, "OIII[%d] = % 0.16g\n",
- i, OIII);*/
- j++;
-
- OIV = input[j];
- /*fprintf(stderr, "OIV[%d] = % 0.16g\n",
- i, OIV);*/
- j++;
-
- OV = input[j];
- /*fprintf(stderr, "OV[%d] = % 0.16g\n",
- i, OV);*/
- j++;
-
- OVI = input[j];
- /*fprintf(stderr, "OVI[%d] = % 0.16g\n",
- i, OVI);*/
- j++;
-
- OVII = input[j];
- /*fprintf(stderr, "OVII[%d] = % 0.16g\n",
- i, OVII);*/
- j++;
-
- OVIII = input[j];
- /*fprintf(stderr, "OVIII[%d] = % 0.16g\n",
- i, OVIII);*/
- j++;
-
- OIX = input[j];
- /*fprintf(stderr, "OIX[%d] = % 0.16g\n",
- i, OIX);*/
- j++;
-
- density = 1.00794*HI + 1.00794*HII + 4.002602*HeI + 4.002602*HeII + 4.002602*HeIII + 15.9994*OI + 15.9994*OII + 15.9994*OIII + 15.9994*OIV + 15.9994*OIX + 15.9994*OV + 15.9994*OVI + 15.9994*OVII + 15.9994*OVIII;
- data->Ts[i] = density*ge*mh/(kb*(HI/(gamma - 1.0) + HII/(gamma - 1.0) + HeI/(gamma - 1.0) + HeII/(gamma - 1.0) + HeIII/(gamma - 1.0) + OI/(gamma - 1.0) + OII/(gamma - 1.0) + OIII/(gamma - 1.0) + OIV/(gamma - 1.0) + OIX/(gamma - 1.0) + OV/(gamma - 1.0) + OVI/(gamma - 1.0) + OVII/(gamma - 1.0) + OVIII/(gamma - 1.0) + de/(gamma - 1.0)));
- if (data->Ts[i] < data->bounds[0]) {
- data->Ts[i] = data->bounds[0];
- } else if (data->Ts[i] > data->bounds[1]) {
- data->Ts[i] = data->bounds[1];
- }
- data->logTs[i] = log(data->Ts[i]);
- data->invTs[i] = 1.0 / data->Ts[i];
- data->dTs_ge[i] =
- density*mh/(kb*(HI/(gamma - 1.0) + HII/(gamma - 1.0) + HeI/(gamma - 1.0) + HeII/(gamma - 1.0) + HeIII/(gamma - 1.0) + OI/(gamma - 1.0) + OII/(gamma - 1.0) + OIII/(gamma - 1.0) + OIV/(gamma - 1.0) + OIX/(gamma - 1.0) + OV/(gamma - 1.0) + OVI/(gamma - 1.0) + OVII/(gamma - 1.0) + OVIII/(gamma - 1.0) + de/(gamma - 1.0)));
- /*fprintf(stderr, "T[%d] = % 0.16g, density = % 0.16g\n",
- i, data->Ts[i], density);*/
- }
-
- }
-
- /*
- This setup may be different than the user may anticipate, as a result
- of the lockstep timestep we use for a pencil beam through the grid.
- As such, it accepts the number of things to interpolate and makes
- assumptions about the sizes of the rates.
- */
- /* This also requires no templating other than for the solver name...*/
- void po_interpolate_rates(po_data *data,
- int nstrip)
- {
- int i, bin_id, zbin_id;
- double lb, t1, t2;
- double lbz, z1, z2;
- int no_photo = 0;
- lb = log(data->bounds[0]);
- lbz = log(data->z_bounds[0] + 1.0);
- /*fprintf(stderr, "lb = % 0.16g, ub = % 0.16g\n", lb, ub);*/
- for (i = 0; i < nstrip; i++) {
- data->bin_id[i] = bin_id = (int) (data->idbin * (data->logTs[i] - lb));
- if (data->bin_id[i] <= 0) {
- data->bin_id[i] = 0;
- } else if (data->bin_id[i] >= data->nbins) {
- data->bin_id[i] = data->nbins - 1;
- }
- t1 = (lb + (bin_id ) * data->dbin);
- t2 = (lb + (bin_id + 1) * data->dbin);
- data->Tdef[i] = (data->logTs[i] - t1)/(t2 - t1);
- data->dT[i] = (t2 - t1);
- /*fprintf(stderr, "INTERP: %d, bin_id = %d, dT = % 0.16g, T = % 0.16g, logT = % 0.16g\n",
- i, data->bin_id[i], data->dT[i], data->Ts[i],
- data->logTs[i]);*/
- }
-
- if ((data->current_z >= data->z_bounds[0]) && (data->current_z < data->z_bounds[1])) {
- zbin_id = (int) (data->id_zbin * (log(data->current_z + 1.0) - lbz));
- if (zbin_id <= 0) {
- zbin_id = 0;
- } else if (zbin_id >= data->n_zbins) {
- zbin_id = data->n_zbins - 1;
- }
- z1 = (lbz + (zbin_id ) * data->d_zbin);
- z2 = (lbz + (zbin_id + 1) * data->d_zbin);
- data->zdef = (log(data->current_z + 1.0) - z1)/(z2 - z1);
- data->dz = (exp(z2) - exp(z1)); //note: given this, we don't have to divide rate of change by z
- } else {
- no_photo = 1;
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_HeI_i[i] = data->r_HeI_i[bin_id] +
- data->Tdef[i] * (data->r_HeI_i[bin_id+1] - data->r_HeI_i[bin_id]);
- data->drs_HeI_i[i] = (data->r_HeI_i[bin_id+1] - data->r_HeI_i[bin_id]);
- data->drs_HeI_i[i] /= data->dT[i];
- data->drs_HeI_i[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->rs_HeI_pi[i] = 0.0;
- data->drs_HeI_pi[i] = 0.0;
- } else {
- data->rs_HeI_pi[i] = data->r_HeI_pi[zbin_id] +
- data->zdef * (data->r_HeI_pi[zbin_id+1] - data->r_HeI_pi[zbin_id]);
- data->drs_HeI_pi[i] = (data->r_HeI_pi[zbin_id+1] - data->r_HeI_pi[zbin_id]);
- data->drs_HeI_pi[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_HeII_i[i] = data->r_HeII_i[bin_id] +
- data->Tdef[i] * (data->r_HeII_i[bin_id+1] - data->r_HeII_i[bin_id]);
- data->drs_HeII_i[i] = (data->r_HeII_i[bin_id+1] - data->r_HeII_i[bin_id]);
- data->drs_HeII_i[i] /= data->dT[i];
- data->drs_HeII_i[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->rs_HeII_pi[i] = 0.0;
- data->drs_HeII_pi[i] = 0.0;
- } else {
- data->rs_HeII_pi[i] = data->r_HeII_pi[zbin_id] +
- data->zdef * (data->r_HeII_pi[zbin_id+1] - data->r_HeII_pi[zbin_id]);
- data->drs_HeII_pi[i] = (data->r_HeII_pi[zbin_id+1] - data->r_HeII_pi[zbin_id]);
- data->drs_HeII_pi[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_HeII_r[i] = data->r_HeII_r[bin_id] +
- data->Tdef[i] * (data->r_HeII_r[bin_id+1] - data->r_HeII_r[bin_id]);
- data->drs_HeII_r[i] = (data->r_HeII_r[bin_id+1] - data->r_HeII_r[bin_id]);
- data->drs_HeII_r[i] /= data->dT[i];
- data->drs_HeII_r[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_HeIII_r[i] = data->r_HeIII_r[bin_id] +
- data->Tdef[i] * (data->r_HeIII_r[bin_id+1] - data->r_HeIII_r[bin_id]);
- data->drs_HeIII_r[i] = (data->r_HeIII_r[bin_id+1] - data->r_HeIII_r[bin_id]);
- data->drs_HeIII_r[i] /= data->dT[i];
- data->drs_HeIII_r[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_HI_i[i] = data->r_HI_i[bin_id] +
- data->Tdef[i] * (data->r_HI_i[bin_id+1] - data->r_HI_i[bin_id]);
- data->drs_HI_i[i] = (data->r_HI_i[bin_id+1] - data->r_HI_i[bin_id]);
- data->drs_HI_i[i] /= data->dT[i];
- data->drs_HI_i[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->rs_HI_pi[i] = 0.0;
- data->drs_HI_pi[i] = 0.0;
- } else {
- data->rs_HI_pi[i] = data->r_HI_pi[zbin_id] +
- data->zdef * (data->r_HI_pi[zbin_id+1] - data->r_HI_pi[zbin_id]);
- data->drs_HI_pi[i] = (data->r_HI_pi[zbin_id+1] - data->r_HI_pi[zbin_id]);
- data->drs_HI_pi[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_HII_r[i] = data->r_HII_r[bin_id] +
- data->Tdef[i] * (data->r_HII_r[bin_id+1] - data->r_HII_r[bin_id]);
- data->drs_HII_r[i] = (data->r_HII_r[bin_id+1] - data->r_HII_r[bin_id]);
- data->drs_HII_r[i] /= data->dT[i];
- data->drs_HII_r[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OI_i[i] = data->r_OI_i[bin_id] +
- data->Tdef[i] * (data->r_OI_i[bin_id+1] - data->r_OI_i[bin_id]);
- data->drs_OI_i[i] = (data->r_OI_i[bin_id+1] - data->r_OI_i[bin_id]);
- data->drs_OI_i[i] /= data->dT[i];
- data->drs_OI_i[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->rs_OI_pi[i] = 0.0;
- data->drs_OI_pi[i] = 0.0;
- } else {
- data->rs_OI_pi[i] = data->r_OI_pi[zbin_id] +
- data->zdef * (data->r_OI_pi[zbin_id+1] - data->r_OI_pi[zbin_id]);
- data->drs_OI_pi[i] = (data->r_OI_pi[zbin_id+1] - data->r_OI_pi[zbin_id]);
- data->drs_OI_pi[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OII_i[i] = data->r_OII_i[bin_id] +
- data->Tdef[i] * (data->r_OII_i[bin_id+1] - data->r_OII_i[bin_id]);
- data->drs_OII_i[i] = (data->r_OII_i[bin_id+1] - data->r_OII_i[bin_id]);
- data->drs_OII_i[i] /= data->dT[i];
- data->drs_OII_i[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->rs_OII_pi[i] = 0.0;
- data->drs_OII_pi[i] = 0.0;
- } else {
- data->rs_OII_pi[i] = data->r_OII_pi[zbin_id] +
- data->zdef * (data->r_OII_pi[zbin_id+1] - data->r_OII_pi[zbin_id]);
- data->drs_OII_pi[i] = (data->r_OII_pi[zbin_id+1] - data->r_OII_pi[zbin_id]);
- data->drs_OII_pi[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OII_r[i] = data->r_OII_r[bin_id] +
- data->Tdef[i] * (data->r_OII_r[bin_id+1] - data->r_OII_r[bin_id]);
- data->drs_OII_r[i] = (data->r_OII_r[bin_id+1] - data->r_OII_r[bin_id]);
- data->drs_OII_r[i] /= data->dT[i];
- data->drs_OII_r[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OIII_i[i] = data->r_OIII_i[bin_id] +
- data->Tdef[i] * (data->r_OIII_i[bin_id+1] - data->r_OIII_i[bin_id]);
- data->drs_OIII_i[i] = (data->r_OIII_i[bin_id+1] - data->r_OIII_i[bin_id]);
- data->drs_OIII_i[i] /= data->dT[i];
- data->drs_OIII_i[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->rs_OIII_pi[i] = 0.0;
- data->drs_OIII_pi[i] = 0.0;
- } else {
- data->rs_OIII_pi[i] = data->r_OIII_pi[zbin_id] +
- data->zdef * (data->r_OIII_pi[zbin_id+1] - data->r_OIII_pi[zbin_id]);
- data->drs_OIII_pi[i] = (data->r_OIII_pi[zbin_id+1] - data->r_OIII_pi[zbin_id]);
- data->drs_OIII_pi[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OIII_r[i] = data->r_OIII_r[bin_id] +
- data->Tdef[i] * (data->r_OIII_r[bin_id+1] - data->r_OIII_r[bin_id]);
- data->drs_OIII_r[i] = (data->r_OIII_r[bin_id+1] - data->r_OIII_r[bin_id]);
- data->drs_OIII_r[i] /= data->dT[i];
- data->drs_OIII_r[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OIV_i[i] = data->r_OIV_i[bin_id] +
- data->Tdef[i] * (data->r_OIV_i[bin_id+1] - data->r_OIV_i[bin_id]);
- data->drs_OIV_i[i] = (data->r_OIV_i[bin_id+1] - data->r_OIV_i[bin_id]);
- data->drs_OIV_i[i] /= data->dT[i];
- data->drs_OIV_i[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->rs_OIV_pi[i] = 0.0;
- data->drs_OIV_pi[i] = 0.0;
- } else {
- data->rs_OIV_pi[i] = data->r_OIV_pi[zbin_id] +
- data->zdef * (data->r_OIV_pi[zbin_id+1] - data->r_OIV_pi[zbin_id]);
- data->drs_OIV_pi[i] = (data->r_OIV_pi[zbin_id+1] - data->r_OIV_pi[zbin_id]);
- data->drs_OIV_pi[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OIV_r[i] = data->r_OIV_r[bin_id] +
- data->Tdef[i] * (data->r_OIV_r[bin_id+1] - data->r_OIV_r[bin_id]);
- data->drs_OIV_r[i] = (data->r_OIV_r[bin_id+1] - data->r_OIV_r[bin_id]);
- data->drs_OIV_r[i] /= data->dT[i];
- data->drs_OIV_r[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OIX_r[i] = data->r_OIX_r[bin_id] +
- data->Tdef[i] * (data->r_OIX_r[bin_id+1] - data->r_OIX_r[bin_id]);
- data->drs_OIX_r[i] = (data->r_OIX_r[bin_id+1] - data->r_OIX_r[bin_id]);
- data->drs_OIX_r[i] /= data->dT[i];
- data->drs_OIX_r[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OV_i[i] = data->r_OV_i[bin_id] +
- data->Tdef[i] * (data->r_OV_i[bin_id+1] - data->r_OV_i[bin_id]);
- data->drs_OV_i[i] = (data->r_OV_i[bin_id+1] - data->r_OV_i[bin_id]);
- data->drs_OV_i[i] /= data->dT[i];
- data->drs_OV_i[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->rs_OV_pi[i] = 0.0;
- data->drs_OV_pi[i] = 0.0;
- } else {
- data->rs_OV_pi[i] = data->r_OV_pi[zbin_id] +
- data->zdef * (data->r_OV_pi[zbin_id+1] - data->r_OV_pi[zbin_id]);
- data->drs_OV_pi[i] = (data->r_OV_pi[zbin_id+1] - data->r_OV_pi[zbin_id]);
- data->drs_OV_pi[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OV_r[i] = data->r_OV_r[bin_id] +
- data->Tdef[i] * (data->r_OV_r[bin_id+1] - data->r_OV_r[bin_id]);
- data->drs_OV_r[i] = (data->r_OV_r[bin_id+1] - data->r_OV_r[bin_id]);
- data->drs_OV_r[i] /= data->dT[i];
- data->drs_OV_r[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OVI_i[i] = data->r_OVI_i[bin_id] +
- data->Tdef[i] * (data->r_OVI_i[bin_id+1] - data->r_OVI_i[bin_id]);
- data->drs_OVI_i[i] = (data->r_OVI_i[bin_id+1] - data->r_OVI_i[bin_id]);
- data->drs_OVI_i[i] /= data->dT[i];
- data->drs_OVI_i[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->rs_OVI_pi[i] = 0.0;
- data->drs_OVI_pi[i] = 0.0;
- } else {
- data->rs_OVI_pi[i] = data->r_OVI_pi[zbin_id] +
- data->zdef * (data->r_OVI_pi[zbin_id+1] - data->r_OVI_pi[zbin_id]);
- data->drs_OVI_pi[i] = (data->r_OVI_pi[zbin_id+1] - data->r_OVI_pi[zbin_id]);
- data->drs_OVI_pi[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OVI_r[i] = data->r_OVI_r[bin_id] +
- data->Tdef[i] * (data->r_OVI_r[bin_id+1] - data->r_OVI_r[bin_id]);
- data->drs_OVI_r[i] = (data->r_OVI_r[bin_id+1] - data->r_OVI_r[bin_id]);
- data->drs_OVI_r[i] /= data->dT[i];
- data->drs_OVI_r[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OVII_i[i] = data->r_OVII_i[bin_id] +
- data->Tdef[i] * (data->r_OVII_i[bin_id+1] - data->r_OVII_i[bin_id]);
- data->drs_OVII_i[i] = (data->r_OVII_i[bin_id+1] - data->r_OVII_i[bin_id]);
- data->drs_OVII_i[i] /= data->dT[i];
- data->drs_OVII_i[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->rs_OVII_pi[i] = 0.0;
- data->drs_OVII_pi[i] = 0.0;
- } else {
- data->rs_OVII_pi[i] = data->r_OVII_pi[zbin_id] +
- data->zdef * (data->r_OVII_pi[zbin_id+1] - data->r_OVII_pi[zbin_id]);
- data->drs_OVII_pi[i] = (data->r_OVII_pi[zbin_id+1] - data->r_OVII_pi[zbin_id]);
- data->drs_OVII_pi[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OVII_r[i] = data->r_OVII_r[bin_id] +
- data->Tdef[i] * (data->r_OVII_r[bin_id+1] - data->r_OVII_r[bin_id]);
- data->drs_OVII_r[i] = (data->r_OVII_r[bin_id+1] - data->r_OVII_r[bin_id]);
- data->drs_OVII_r[i] /= data->dT[i];
- data->drs_OVII_r[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OVIII_i[i] = data->r_OVIII_i[bin_id] +
- data->Tdef[i] * (data->r_OVIII_i[bin_id+1] - data->r_OVIII_i[bin_id]);
- data->drs_OVIII_i[i] = (data->r_OVIII_i[bin_id+1] - data->r_OVIII_i[bin_id]);
- data->drs_OVIII_i[i] /= data->dT[i];
- data->drs_OVIII_i[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->rs_OVIII_pi[i] = 0.0;
- data->drs_OVIII_pi[i] = 0.0;
- } else {
- data->rs_OVIII_pi[i] = data->r_OVIII_pi[zbin_id] +
- data->zdef * (data->r_OVIII_pi[zbin_id+1] - data->r_OVIII_pi[zbin_id]);
- data->drs_OVIII_pi[i] = (data->r_OVIII_pi[zbin_id+1] - data->r_OVIII_pi[zbin_id]);
- data->drs_OVIII_pi[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->rs_OVIII_r[i] = data->r_OVIII_r[bin_id] +
- data->Tdef[i] * (data->r_OVIII_r[bin_id+1] - data->r_OVIII_r[bin_id]);
- data->drs_OVIII_r[i] = (data->r_OVIII_r[bin_id+1] - data->r_OVIII_r[bin_id]);
- data->drs_OVIII_r[i] /= data->dT[i];
- data->drs_OVIII_r[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_compton_comp[i] = data->c_compton_comp[bin_id] +
- data->Tdef[i] * (data->c_compton_comp[bin_id+1] - data->c_compton_comp[bin_id]);
- data->dcs_compton_comp[i] = (data->c_compton_comp[bin_id+1] - data->c_compton_comp[bin_id]);;
- data->dcs_compton_comp[i] /= data->dT[i];
- data->dcs_compton_comp[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_HeI_c_HeI_c[i] = data->c_HeI_c_HeI_c[bin_id] +
- data->Tdef[i] * (data->c_HeI_c_HeI_c[bin_id+1] - data->c_HeI_c_HeI_c[bin_id]);
- data->dcs_HeI_c_HeI_c[i] = (data->c_HeI_c_HeI_c[bin_id+1] - data->c_HeI_c_HeI_c[bin_id]);;
- data->dcs_HeI_c_HeI_c[i] /= data->dT[i];
- data->dcs_HeI_c_HeI_c[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->cs_HeI_ph_HeI_ph[i] = 0.0;
- data->dcs_HeI_ph_HeI_ph[i] = 0.0;
- } else {
- data->cs_HeI_ph_HeI_ph[i] = data->c_HeI_ph_HeI_ph[zbin_id] +
- data->zdef * (data->c_HeI_ph_HeI_ph[zbin_id+1] - data->c_HeI_ph_HeI_ph[zbin_id]);
- data->dcs_HeI_ph_HeI_ph[i] = (data->c_HeI_ph_HeI_ph[zbin_id+1] - data->c_HeI_ph_HeI_ph[zbin_id]);;
- data->dcs_HeI_ph_HeI_ph[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_HeII_c_HeII_c[i] = data->c_HeII_c_HeII_c[bin_id] +
- data->Tdef[i] * (data->c_HeII_c_HeII_c[bin_id+1] - data->c_HeII_c_HeII_c[bin_id]);
- data->dcs_HeII_c_HeII_c[i] = (data->c_HeII_c_HeII_c[bin_id+1] - data->c_HeII_c_HeII_c[bin_id]);;
- data->dcs_HeII_c_HeII_c[i] /= data->dT[i];
- data->dcs_HeII_c_HeII_c[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->cs_HeII_ph_HeII_ph[i] = 0.0;
- data->dcs_HeII_ph_HeII_ph[i] = 0.0;
- } else {
- data->cs_HeII_ph_HeII_ph[i] = data->c_HeII_ph_HeII_ph[zbin_id] +
- data->zdef * (data->c_HeII_ph_HeII_ph[zbin_id+1] - data->c_HeII_ph_HeII_ph[zbin_id]);
- data->dcs_HeII_ph_HeII_ph[i] = (data->c_HeII_ph_HeII_ph[zbin_id+1] - data->c_HeII_ph_HeII_ph[zbin_id]);;
- data->dcs_HeII_ph_HeII_ph[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_HeIII_c_HeIII_c[i] = data->c_HeIII_c_HeIII_c[bin_id] +
- data->Tdef[i] * (data->c_HeIII_c_HeIII_c[bin_id+1] - data->c_HeIII_c_HeIII_c[bin_id]);
- data->dcs_HeIII_c_HeIII_c[i] = (data->c_HeIII_c_HeIII_c[bin_id+1] - data->c_HeIII_c_HeIII_c[bin_id]);;
- data->dcs_HeIII_c_HeIII_c[i] /= data->dT[i];
- data->dcs_HeIII_c_HeIII_c[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_HI_c_HI_c[i] = data->c_HI_c_HI_c[bin_id] +
- data->Tdef[i] * (data->c_HI_c_HI_c[bin_id+1] - data->c_HI_c_HI_c[bin_id]);
- data->dcs_HI_c_HI_c[i] = (data->c_HI_c_HI_c[bin_id+1] - data->c_HI_c_HI_c[bin_id]);;
- data->dcs_HI_c_HI_c[i] /= data->dT[i];
- data->dcs_HI_c_HI_c[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->cs_HI_ph_HI_ph[i] = 0.0;
- data->dcs_HI_ph_HI_ph[i] = 0.0;
- } else {
- data->cs_HI_ph_HI_ph[i] = data->c_HI_ph_HI_ph[zbin_id] +
- data->zdef * (data->c_HI_ph_HI_ph[zbin_id+1] - data->c_HI_ph_HI_ph[zbin_id]);
- data->dcs_HI_ph_HI_ph[i] = (data->c_HI_ph_HI_ph[zbin_id+1] - data->c_HI_ph_HI_ph[zbin_id]);;
- data->dcs_HI_ph_HI_ph[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_HII_c_HII_c[i] = data->c_HII_c_HII_c[bin_id] +
- data->Tdef[i] * (data->c_HII_c_HII_c[bin_id+1] - data->c_HII_c_HII_c[bin_id]);
- data->dcs_HII_c_HII_c[i] = (data->c_HII_c_HII_c[bin_id+1] - data->c_HII_c_HII_c[bin_id]);;
- data->dcs_HII_c_HII_c[i] /= data->dT[i];
- data->dcs_HII_c_HII_c[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_OI_c_OI_c[i] = data->c_OI_c_OI_c[bin_id] +
- data->Tdef[i] * (data->c_OI_c_OI_c[bin_id+1] - data->c_OI_c_OI_c[bin_id]);
- data->dcs_OI_c_OI_c[i] = (data->c_OI_c_OI_c[bin_id+1] - data->c_OI_c_OI_c[bin_id]);;
- data->dcs_OI_c_OI_c[i] /= data->dT[i];
- data->dcs_OI_c_OI_c[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->cs_OI_ph_OI_ph[i] = 0.0;
- data->dcs_OI_ph_OI_ph[i] = 0.0;
- } else {
- data->cs_OI_ph_OI_ph[i] = data->c_OI_ph_OI_ph[zbin_id] +
- data->zdef * (data->c_OI_ph_OI_ph[zbin_id+1] - data->c_OI_ph_OI_ph[zbin_id]);
- data->dcs_OI_ph_OI_ph[i] = (data->c_OI_ph_OI_ph[zbin_id+1] - data->c_OI_ph_OI_ph[zbin_id]);;
- data->dcs_OI_ph_OI_ph[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_OII_c_OII_c[i] = data->c_OII_c_OII_c[bin_id] +
- data->Tdef[i] * (data->c_OII_c_OII_c[bin_id+1] - data->c_OII_c_OII_c[bin_id]);
- data->dcs_OII_c_OII_c[i] = (data->c_OII_c_OII_c[bin_id+1] - data->c_OII_c_OII_c[bin_id]);;
- data->dcs_OII_c_OII_c[i] /= data->dT[i];
- data->dcs_OII_c_OII_c[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->cs_OII_ph_OII_ph[i] = 0.0;
- data->dcs_OII_ph_OII_ph[i] = 0.0;
- } else {
- data->cs_OII_ph_OII_ph[i] = data->c_OII_ph_OII_ph[zbin_id] +
- data->zdef * (data->c_OII_ph_OII_ph[zbin_id+1] - data->c_OII_ph_OII_ph[zbin_id]);
- data->dcs_OII_ph_OII_ph[i] = (data->c_OII_ph_OII_ph[zbin_id+1] - data->c_OII_ph_OII_ph[zbin_id]);;
- data->dcs_OII_ph_OII_ph[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_OIII_c_OIII_c[i] = data->c_OIII_c_OIII_c[bin_id] +
- data->Tdef[i] * (data->c_OIII_c_OIII_c[bin_id+1] - data->c_OIII_c_OIII_c[bin_id]);
- data->dcs_OIII_c_OIII_c[i] = (data->c_OIII_c_OIII_c[bin_id+1] - data->c_OIII_c_OIII_c[bin_id]);;
- data->dcs_OIII_c_OIII_c[i] /= data->dT[i];
- data->dcs_OIII_c_OIII_c[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->cs_OIII_ph_OIII_ph[i] = 0.0;
- data->dcs_OIII_ph_OIII_ph[i] = 0.0;
- } else {
- data->cs_OIII_ph_OIII_ph[i] = data->c_OIII_ph_OIII_ph[zbin_id] +
- data->zdef * (data->c_OIII_ph_OIII_ph[zbin_id+1] - data->c_OIII_ph_OIII_ph[zbin_id]);
- data->dcs_OIII_ph_OIII_ph[i] = (data->c_OIII_ph_OIII_ph[zbin_id+1] - data->c_OIII_ph_OIII_ph[zbin_id]);;
- data->dcs_OIII_ph_OIII_ph[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_OIV_c_OIV_c[i] = data->c_OIV_c_OIV_c[bin_id] +
- data->Tdef[i] * (data->c_OIV_c_OIV_c[bin_id+1] - data->c_OIV_c_OIV_c[bin_id]);
- data->dcs_OIV_c_OIV_c[i] = (data->c_OIV_c_OIV_c[bin_id+1] - data->c_OIV_c_OIV_c[bin_id]);;
- data->dcs_OIV_c_OIV_c[i] /= data->dT[i];
- data->dcs_OIV_c_OIV_c[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- if (no_photo) {
- data->cs_OIV_ph_OIV_ph[i] = 0.0;
- data->dcs_OIV_ph_OIV_ph[i] = 0.0;
- } else {
- data->cs_OIV_ph_OIV_ph[i] = data->c_OIV_ph_OIV_ph[zbin_id] +
- data->zdef * (data->c_OIV_ph_OIV_ph[zbin_id+1] - data->c_OIV_ph_OIV_ph[zbin_id]);
- data->dcs_OIV_ph_OIV_ph[i] = (data->c_OIV_ph_OIV_ph[zbin_id+1] - data->c_OIV_ph_OIV_ph[zbin_id]);;
- data->dcs_OIV_ph_OIV_ph[i] /= data->dz;
- }
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_OIX_c_OIX_c[i] = data->c_OIX_c_OIX_c[bin_id] +
- data->Tdef[i] * (data->c_OIX_c_OIX_c[bin_id+1] - data->c_OIX_c_OIX_c[bin_id]);
- data->dcs_OIX_c_OIX_c[i] = (data->c_OIX_c_OIX_c[bin_id+1] - data->c_OIX_c_OIX_c[bin_id]);;
- data->dcs_OIX_c_OIX_c[i] /= data->dT[i];
- data->dcs_OIX_c_OIX_c[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- bin_id = data->bin_id[i];
- data->cs_OV_c_OV_c[i] = data->c_OV_c_OV_c[bin_id] +
- data->Tdef[i] * (data->c_OV_c_OV_c[bin_id+1] - data->c_OV_c_OV_c[bin_id]);
- data->dcs_OV_c_OV_c[i] = (data->c_OV_c_OV_c[bin_id+1] - data->c_OV_c_OV_c[bin_id]);;
- data->dcs_OV_c_OV_c[i] /= data->dT[i];
- data->dcs_OV_c_OV_c[i] *= data->invTs[i];
- }
-
- for (i = 0; i < nstrip; i++) {
- …
Large files files are truncated, but you can click here to view the full file