PageRenderTime 407ms CodeModel.GetById 170ms app.highlight 216ms RepoModel.GetById 1ms app.codeStats 3ms

/src/enzo/oxygen_solver.C

https://bitbucket.org/devinsilvia/enzo-dengo-devin
C++ | 2123 lines | 1231 code | 563 blank | 329 comment | 149 complexity | b918f9adf7d459fef04598a0d7a5b104 MD5 | raw file

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

   1
   2/* THIS FILE HAS BEEN AUTO-GENERATED.  DO NOT EDIT. */
   3
   4/* This is C++ code to read HDF5 files for
   5   reaction rates, cooling rates, and initial
   6   conditions for the chemical network defined
   7   by the user.  In addition, this contains
   8   code for calculating temperature from the
   9   gas energy and computing the RHS and the
  10   Jacobian of the system of equations which
  11   will be fed into the solver.
  12*/
  13
  14
  15#include <alloca.h>
  16#include "oxygen_solver.h"
  17
  18oxygen_data *oxygen_setup_data(void) {
  19
  20    oxygen_data *data = (oxygen_data *) malloc(sizeof(oxygen_data));
  21
  22    data->bounds[0] = 1.0;
  23    data->bounds[1] = 1000000000000.0;
  24    data->nbins = 1024;
  25    data->dbin = (log(data->bounds[1]) - log(data->bounds[0])) / data->nbins;
  26    data->idbin = 1.0L / data->dbin;
  27    
  28    oxygen_read_rate_tables(data);
  29    fprintf(stderr, "Successfully read in rate tables.\n");
  30
  31    oxygen_read_cooling_tables(data);
  32    fprintf(stderr, "Successfully read in cooling rate tables.\n");
  33
  34    return data;
  35
  36}
  37
  38
  39int oxygen_main(int argc, char** argv)
  40{
  41    oxygen_data *data = oxygen_setup_data();
  42
  43    /* Initial conditions */
  44
  45    hid_t file_id = H5Fopen("oxygen_initial_conditions.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
  46    if (file_id < 0) {fprintf(stderr, "Failed to open "
  47        "oxygen_initial_conditions.h5 so dying.\n");
  48        return(1);}
  49
  50    /* Allocate the correct number of cells */
  51    hsize_t dims; /* We have flat versus number of species */
  52
  53    /* Check gas energy to get the number of cells */
  54    fprintf(stderr, "Getting dimensionality from ge:\n");
  55    herr_t status = H5LTget_dataset_info(file_id, "/ge", &dims, NULL, NULL);
  56    if(status == -1) {
  57        fprintf(stderr, "Error opening initial conditions file.\n");
  58        return 1;
  59    }
  60    fprintf(stderr, "  ncells = % 3i\n", (int) dims);
  61    data->ncells = dims;
  62
  63    int N = 11;
  64
  65    double *atol, *rtol;
  66    atol = (double *) alloca(N * dims * sizeof(double));
  67    rtol = (double *) alloca(N * dims * sizeof(double));
  68
  69    double *tics = (double *) alloca(dims * sizeof(double));
  70    double *ics = (double *) alloca(dims * N * sizeof(double));
  71    double *input = (double *) alloca(dims * N * sizeof(double));
  72    
  73    unsigned int i = 0, j;
  74    
  75    fprintf(stderr, "Reading I.C. for /ge\n");
  76    H5LTread_dataset_double(file_id, "/ge", tics);
  77    for (j = 0; j < dims; j++) {
  78        ics[j * N + i] = tics[j] / 1.0; /* Convert to number density */
  79        atol[j * N + i] = tics[j] * 1e-06;
  80        rtol[j * N + i] = 1e-06;
  81        if(j==0) {
  82            fprintf(stderr, "ge[0] = %0.3g, atol => % 0.16g\n",
  83                    tics[j], atol[j]);
  84        }
  85    }
  86    i++;
  87    
  88    fprintf(stderr, "Reading I.C. for /de\n");
  89    H5LTread_dataset_double(file_id, "/de", tics);
  90    for (j = 0; j < dims; j++) {
  91        ics[j * N + i] = tics[j] / 1.0; /* Convert to number density */
  92        atol[j * N + i] = tics[j] * 1e-06;
  93        rtol[j * N + i] = 1e-06;
  94        if(j==0) {
  95            fprintf(stderr, "de[0] = %0.3g, atol => % 0.16g\n",
  96                    tics[j], atol[j]);
  97        }
  98    }
  99    i++;
 100    
 101    fprintf(stderr, "Reading I.C. for /OI\n");
 102    H5LTread_dataset_double(file_id, "/OI", tics);
 103    for (j = 0; j < dims; j++) {
 104        ics[j * N + i] = tics[j] / 16; /* Convert to number density */
 105        atol[j * N + i] = tics[j] * 1e-06;
 106        rtol[j * N + i] = 1e-06;
 107        if(j==0) {
 108            fprintf(stderr, "OI[0] = %0.3g, atol => % 0.16g\n",
 109                    tics[j], atol[j]);
 110        }
 111    }
 112    i++;
 113    
 114    fprintf(stderr, "Reading I.C. for /OVII\n");
 115    H5LTread_dataset_double(file_id, "/OVII", tics);
 116    for (j = 0; j < dims; j++) {
 117        ics[j * N + i] = tics[j] / 16; /* Convert to number density */
 118        atol[j * N + i] = tics[j] * 1e-06;
 119        rtol[j * N + i] = 1e-06;
 120        if(j==0) {
 121            fprintf(stderr, "OVII[0] = %0.3g, atol => % 0.16g\n",
 122                    tics[j], atol[j]);
 123        }
 124    }
 125    i++;
 126    
 127    fprintf(stderr, "Reading I.C. for /OII\n");
 128    H5LTread_dataset_double(file_id, "/OII", tics);
 129    for (j = 0; j < dims; j++) {
 130        ics[j * N + i] = tics[j] / 16; /* Convert to number density */
 131        atol[j * N + i] = tics[j] * 1e-06;
 132        rtol[j * N + i] = 1e-06;
 133        if(j==0) {
 134            fprintf(stderr, "OII[0] = %0.3g, atol => % 0.16g\n",
 135                    tics[j], atol[j]);
 136        }
 137    }
 138    i++;
 139    
 140    fprintf(stderr, "Reading I.C. for /OVI\n");
 141    H5LTread_dataset_double(file_id, "/OVI", tics);
 142    for (j = 0; j < dims; j++) {
 143        ics[j * N + i] = tics[j] / 16; /* Convert to number density */
 144        atol[j * N + i] = tics[j] * 1e-06;
 145        rtol[j * N + i] = 1e-06;
 146        if(j==0) {
 147            fprintf(stderr, "OVI[0] = %0.3g, atol => % 0.16g\n",
 148                    tics[j], atol[j]);
 149        }
 150    }
 151    i++;
 152    
 153    fprintf(stderr, "Reading I.C. for /OIII\n");
 154    H5LTread_dataset_double(file_id, "/OIII", tics);
 155    for (j = 0; j < dims; j++) {
 156        ics[j * N + i] = tics[j] / 16; /* Convert to number density */
 157        atol[j * N + i] = tics[j] * 1e-06;
 158        rtol[j * N + i] = 1e-06;
 159        if(j==0) {
 160            fprintf(stderr, "OIII[0] = %0.3g, atol => % 0.16g\n",
 161                    tics[j], atol[j]);
 162        }
 163    }
 164    i++;
 165    
 166    fprintf(stderr, "Reading I.C. for /OIV\n");
 167    H5LTread_dataset_double(file_id, "/OIV", tics);
 168    for (j = 0; j < dims; j++) {
 169        ics[j * N + i] = tics[j] / 16; /* Convert to number density */
 170        atol[j * N + i] = tics[j] * 1e-06;
 171        rtol[j * N + i] = 1e-06;
 172        if(j==0) {
 173            fprintf(stderr, "OIV[0] = %0.3g, atol => % 0.16g\n",
 174                    tics[j], atol[j]);
 175        }
 176    }
 177    i++;
 178    
 179    fprintf(stderr, "Reading I.C. for /OV\n");
 180    H5LTread_dataset_double(file_id, "/OV", tics);
 181    for (j = 0; j < dims; j++) {
 182        ics[j * N + i] = tics[j] / 16; /* Convert to number density */
 183        atol[j * N + i] = tics[j] * 1e-06;
 184        rtol[j * N + i] = 1e-06;
 185        if(j==0) {
 186            fprintf(stderr, "OV[0] = %0.3g, atol => % 0.16g\n",
 187                    tics[j], atol[j]);
 188        }
 189    }
 190    i++;
 191    
 192    fprintf(stderr, "Reading I.C. for /OVIII\n");
 193    H5LTread_dataset_double(file_id, "/OVIII", tics);
 194    for (j = 0; j < dims; j++) {
 195        ics[j * N + i] = tics[j] / 16; /* Convert to number density */
 196        atol[j * N + i] = tics[j] * 1e-06;
 197        rtol[j * N + i] = 1e-06;
 198        if(j==0) {
 199            fprintf(stderr, "OVIII[0] = %0.3g, atol => % 0.16g\n",
 200                    tics[j], atol[j]);
 201        }
 202    }
 203    i++;
 204    
 205    fprintf(stderr, "Reading I.C. for /OIX\n");
 206    H5LTread_dataset_double(file_id, "/OIX", tics);
 207    for (j = 0; j < dims; j++) {
 208        ics[j * N + i] = tics[j] / 16; /* Convert to number density */
 209        atol[j * N + i] = tics[j] * 1e-06;
 210        rtol[j * N + i] = 1e-06;
 211        if(j==0) {
 212            fprintf(stderr, "OIX[0] = %0.3g, atol => % 0.16g\n",
 213                    tics[j], atol[j]);
 214        }
 215    }
 216    i++;
 217    
 218
 219    H5Fclose(file_id);
 220
 221    double dtf = 2e17;
 222    double dt = -1.0;
 223    for (i = 0; i < dims * N; i++) input[i] = ics[i];
 224    double ttot;
 225    ttot = dengo_evolve_oxygen(dtf, dt, input, rtol, atol, dims, data);
 226
 227    /* Write results to HDF5 file */
 228    file_id = H5Fcreate("oxygen_solution.h5", H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
 229    hsize_t dimsarr[1];
 230    dimsarr[0] = dims;
 231    i = 0;
 232    
 233    double ge[dims];
 234    for (j = 0; j < dims; j++) {
 235        ge[j] = input[j * N + i] * 1.0;
 236    }
 237    fprintf(stderr, "Writing solution for /ge\n");
 238    H5LTmake_dataset_double(file_id, "/ge", 1, dimsarr, ge);
 239    i++;
 240    
 241    double de[dims];
 242    for (j = 0; j < dims; j++) {
 243        de[j] = input[j * N + i] * 1.0;
 244    }
 245    fprintf(stderr, "Writing solution for /de\n");
 246    H5LTmake_dataset_double(file_id, "/de", 1, dimsarr, de);
 247    i++;
 248    
 249    double OI[dims];
 250    for (j = 0; j < dims; j++) {
 251        OI[j] = input[j * N + i] * 16;
 252    }
 253    fprintf(stderr, "Writing solution for /OI\n");
 254    H5LTmake_dataset_double(file_id, "/OI", 1, dimsarr, OI);
 255    i++;
 256    
 257    double OVII[dims];
 258    for (j = 0; j < dims; j++) {
 259        OVII[j] = input[j * N + i] * 16;
 260    }
 261    fprintf(stderr, "Writing solution for /OVII\n");
 262    H5LTmake_dataset_double(file_id, "/OVII", 1, dimsarr, OVII);
 263    i++;
 264    
 265    double OII[dims];
 266    for (j = 0; j < dims; j++) {
 267        OII[j] = input[j * N + i] * 16;
 268    }
 269    fprintf(stderr, "Writing solution for /OII\n");
 270    H5LTmake_dataset_double(file_id, "/OII", 1, dimsarr, OII);
 271    i++;
 272    
 273    double OVI[dims];
 274    for (j = 0; j < dims; j++) {
 275        OVI[j] = input[j * N + i] * 16;
 276    }
 277    fprintf(stderr, "Writing solution for /OVI\n");
 278    H5LTmake_dataset_double(file_id, "/OVI", 1, dimsarr, OVI);
 279    i++;
 280    
 281    double OIII[dims];
 282    for (j = 0; j < dims; j++) {
 283        OIII[j] = input[j * N + i] * 16;
 284    }
 285    fprintf(stderr, "Writing solution for /OIII\n");
 286    H5LTmake_dataset_double(file_id, "/OIII", 1, dimsarr, OIII);
 287    i++;
 288    
 289    double OIV[dims];
 290    for (j = 0; j < dims; j++) {
 291        OIV[j] = input[j * N + i] * 16;
 292    }
 293    fprintf(stderr, "Writing solution for /OIV\n");
 294    H5LTmake_dataset_double(file_id, "/OIV", 1, dimsarr, OIV);
 295    i++;
 296    
 297    double OV[dims];
 298    for (j = 0; j < dims; j++) {
 299        OV[j] = input[j * N + i] * 16;
 300    }
 301    fprintf(stderr, "Writing solution for /OV\n");
 302    H5LTmake_dataset_double(file_id, "/OV", 1, dimsarr, OV);
 303    i++;
 304    
 305    double OVIII[dims];
 306    for (j = 0; j < dims; j++) {
 307        OVIII[j] = input[j * N + i] * 16;
 308    }
 309    fprintf(stderr, "Writing solution for /OVIII\n");
 310    H5LTmake_dataset_double(file_id, "/OVIII", 1, dimsarr, OVIII);
 311    i++;
 312    
 313    double OIX[dims];
 314    for (j = 0; j < dims; j++) {
 315        OIX[j] = input[j * N + i] * 16;
 316    }
 317    fprintf(stderr, "Writing solution for /OIX\n");
 318    H5LTmake_dataset_double(file_id, "/OIX", 1, dimsarr, OIX);
 319    i++;
 320    
 321    double temperature[dims];
 322    for (j = 0; j < dims; j++) {
 323    	temperature[j] = data->Ts[j];
 324    }
 325    H5LTmake_dataset_double(file_id, "/T", 1, dimsarr, temperature);
 326    double time[1];
 327    time[0] = ttot;
 328    double timestep[1];
 329    timestep[0] = dt;
 330    H5LTset_attribute_double(file_id, "/", "time", time, 1); 
 331    H5LTset_attribute_double(file_id, "/", "timestep", timestep, 1);
 332    H5Fclose(file_id);
 333    
 334    return 0;
 335}
 336 
 337
 338
 339
 340double dengo_evolve_oxygen (double dtf, double &dt, double *input,
 341            double *rtol, double *atol, long long dims, oxygen_data *data) {
 342    int i, j;
 343    hid_t file_id;
 344    /*fprintf(stderr, "  ncells = % 3i\n", (int) dims);*/
 345
 346    int N = 11;
 347    rhs_f f = calculate_rhs_oxygen;
 348    jac_f jf = calculate_jacobian_oxygen;
 349    if (dt < 0) dt = dtf / 1000.0;
 350    int niter = 0;
 351    int siter = 0;
 352    double ttot = 0;
 353    double *scale = (double *) alloca(dims * N * sizeof(double));
 354    double *prev = (double *) alloca(dims * N * sizeof(double));
 355    for (i = 0; i < dims * N; i++) scale[i] = 1.0;
 356    for (i = 0; i < dims * N; i++) prev[i] = input[i];
 357    double *u0 = (double *) alloca(N*dims*sizeof(double));
 358    double *s  = (double *) alloca(N*sizeof(double));
 359    double *gu = (double *) alloca(N*dims*sizeof(double));
 360    double *Ju = (double *) alloca(N*N*dims*sizeof(double));
 361    double floor_value = 1e-25;
 362    double pssum, issum, sumdiff;
 363    int snum;
 364    int relaxed0 = 0;
 365    int relaxed1 = 0;
 366    int relaxed2 = 0;
 367    int maxspeciesi;
 368    double maxspeciesval;
 369
 370    /*
 371    for (i = 0; i < dims * N; i++) {
 372        fprintf(stderr, "rtol[%d] = %0.5g \n", i, rtol[i]);
 373    }
 374    */
 375
 376    while (ttot < dtf) {
 377        int rv = BE_chem_solve(f, jf, input, dt, rtol, atol, dims, N,
 378                               scale, (void *) data, u0, s, gu, Ju);
 379	/*
 380        fprintf(stderr, "Return value [%d]: %i.  %0.5g / %0.5g = %0.5g (%0.5g)\n",
 381                niter, rv, ttot, dtf, ttot/dtf, dt);
 382        fprintf(stderr, "Value[80] = %0.5g %0.5g %0.5g\n",
 383                input[80], prev[80]);//, ics[80]);
 384	*/
 385        for (i = 0; i < dims * N; i++) {
 386            if (input[i] < 0) {
 387                rv = 1;
 388		//fprintf(stderr, "rv = 1, input[%d] = %0.5g\n", i, input[i]);
 389                break;
 390            }
 391        }
 392        if (rv == 0) {
 393            if (siter == 99999999) {
 394	            fprintf(stderr, "Hit maximum number of successful iterations allowed.\n");
 395	            break;
 396            }
 397            siter++;
 398            if (siter % 100000 == 0) {
 399                fprintf(stderr, "Successful Iteration[%d]: (%0.4g) %0.6g / %0.6g = %0.6g \n",
 400                        siter, dt, ttot, dtf, ttot/dtf);
 401            }
 402
 403            ttot += dt;
 404            dt = DMIN(dt * 1.1, dtf - ttot);
 405
 406            if (dt/dtf > 1e-9 && relaxed2 == 1) {
 407                relaxed2 = 0;
 408                for (i = 0; i < dims * N; i++) rtol[i] *= 0.2;
 409            }
 410
 411            if (dt/dtf > 1e-6 && relaxed1 == 1) {
 412                relaxed1 = 0;
 413                for (i = 0; i < dims * N; i++) rtol[i] *= 0.1;
 414            }
 415            
 416            if (dt/dtf > 1e-3 && relaxed0 == 1) {
 417                relaxed0 = 0;
 418                for (i = 0; i < dims * N; i++) rtol[i] *= 0.1;
 419            }
 420
 421            for (i = 0; i < dims * N; i++) prev[i] = input[i];
 422	
 423            for (i = 0; i < dims * N; i++) {
 424                if (input[i] < floor_value)
 425                    input[i] = floor_value; //Do I need to worry about conservation of species
 426                                            //here? dealt with below (maybe)
 427            }
 428
 429            pssum = 0.0;
 430            issum = 0.0;
 431            maxspeciesval = 0.0;
 432            maxspeciesi = -1;
 433            for (i = 0; i < dims * N; i++) { //this conservation check method only works
 434                                             //if solving one cell at a time
 435                snum = i % N;
 436                if (snum == 0) continue;
 437                if (snum == 1) {
 438                    input[i] = prev[i]; //we need to make sure electrons aren't floored
 439                                        //or everything gets weird (I think)
 440                    continue;
 441                }
 442
 443                if (prev[i] > maxspeciesval) { //which species is largest?
 444                  maxspeciesval = prev[i];
 445                  maxspeciesi = i;
 446                }
 447
 448                pssum += prev[i];
 449                issum += input[i];
 450                if (snum == (N - 1)) {
 451                    sumdiff = issum - pssum; //how much extra species density did we
 452                                             //generate by flooring?
 453                    if (sumdiff > 0.0) {
 454                      //fprintf(stderr, "The amount of extra species generated by flooring values was: %0.5g, max species = %d\n", sumdiff, maxspeciesi);
 455                        input[maxspeciesi] -= sumdiff; // to maintain conservation, take it out of
 456                                                           //the most abundant ion (ONLY WORKS IF THERE IS A
 457                                                           //SINGLE ELEMENT)
 458                    }
 459                    pssum = 0.0; // at end, reset
 460                    issum = 0.0; // at end, reset
 461                }
 462            }
 463        } else {
 464
 465            dt /= 2.0;
 466
 467            if (dt/dtf < 1e-6 && relaxed0 == 0) {
 468                relaxed0 = 1;
 469                for (i = 0; i < dims * N; i++) rtol[i] *= 10.0;
 470            }
 471            
 472            if (dt/dtf < 1e-9 && relaxed1 == 0) {
 473                relaxed1 = 1;
 474                for (i = 0; i < dims * N; i++) rtol[i] *= 10.0;
 475            }
 476
 477            if (dt/dtf < 1e-12 && relaxed2 == 0) {
 478                relaxed2 = 1;
 479                for (i = 0; i < dims * N; i++) rtol[i] *= 5.0;
 480            }
 481
 482            for (i = 0; i < dims * N; i++) input[i] = prev[i];
 483            
 484            if (dt/dtf < 1e-15)  {
 485              fprintf(stderr, "Dying! dt/dtf = %0.5g\n", dt/dtf);
 486                break;
 487            }
 488        }
 489        niter++;
 490    }
 491    /*fprintf(stderr, "End: %0.5g / %0.5g (%0.5g)\n",
 492        ttot, dtf, dtf-ttot);*/
 493
 494    return ttot;
 495}
 496 
 497
 498
 499void oxygen_read_rate_tables(oxygen_data *data)
 500{
 501    hid_t file_id = H5Fopen("oxygen_tables.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
 502    /* Allocate the correct number of rate tables */
 503    H5LTread_dataset_double(file_id, "/OI_i", data->r_OI_i);
 504    H5LTread_dataset_double(file_id, "/OII_i", data->r_OII_i);
 505    H5LTread_dataset_double(file_id, "/OII_r", data->r_OII_r);
 506    H5LTread_dataset_double(file_id, "/OIII_i", data->r_OIII_i);
 507    H5LTread_dataset_double(file_id, "/OIII_r", data->r_OIII_r);
 508    H5LTread_dataset_double(file_id, "/OIV_i", data->r_OIV_i);
 509    H5LTread_dataset_double(file_id, "/OIV_r", data->r_OIV_r);
 510    H5LTread_dataset_double(file_id, "/OIX_r", data->r_OIX_r);
 511    H5LTread_dataset_double(file_id, "/OV_i", data->r_OV_i);
 512    H5LTread_dataset_double(file_id, "/OV_r", data->r_OV_r);
 513    H5LTread_dataset_double(file_id, "/OVI_i", data->r_OVI_i);
 514    H5LTread_dataset_double(file_id, "/OVI_r", data->r_OVI_r);
 515    H5LTread_dataset_double(file_id, "/OVII_i", data->r_OVII_i);
 516    H5LTread_dataset_double(file_id, "/OVII_r", data->r_OVII_r);
 517    H5LTread_dataset_double(file_id, "/OVIII_i", data->r_OVIII_i);
 518    H5LTread_dataset_double(file_id, "/OVIII_r", data->r_OVIII_r);
 519
 520    H5Fclose(file_id);
 521}
 522
 523void oxygen_read_cooling_tables(oxygen_data *data)
 524{
 525
 526    hid_t file_id = H5Fopen("oxygen_tables.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
 527    /* Allocate the correct number of rate tables */
 528
 529    H5Fclose(file_id);
 530}
 531
 532 
 533
 534
 535void oxygen_calculate_temperature(oxygen_data *data,
 536                        double *input, int nstrip, int nchem)
 537{
 538    int i, j;
 539    double density;
 540    double kb = 1.3806504e-16; // Boltzmann constant [erg/K]
 541    double mh = 1.67e-24;
 542    double gamma = 5.e0/3.e0;
 543    
 544    /* Calculate total density */
 545    double ge;
 546    double de;
 547    double OI;
 548    double OVII;
 549    double OII;
 550    double OVI;
 551    double OIII;
 552    double OIV;
 553    double OV;
 554    double OVIII;
 555    double OIX;
 556
 557    for (i = 0; i<nstrip; i++) {
 558        j = i * nchem;
 559        ge = input[j];
 560        /*fprintf(stderr, "ge[%d] = % 0.16g\n",
 561                i, ge);*/
 562        j++;
 563    
 564        de = input[j];
 565        /*fprintf(stderr, "de[%d] = % 0.16g\n",
 566                i, de);*/
 567        j++;
 568    
 569        OI = input[j];
 570        /*fprintf(stderr, "OI[%d] = % 0.16g\n",
 571                i, OI);*/
 572        j++;
 573    
 574        OVII = input[j];
 575        /*fprintf(stderr, "OVII[%d] = % 0.16g\n",
 576                i, OVII);*/
 577        j++;
 578    
 579        OII = input[j];
 580        /*fprintf(stderr, "OII[%d] = % 0.16g\n",
 581                i, OII);*/
 582        j++;
 583    
 584        OVI = input[j];
 585        /*fprintf(stderr, "OVI[%d] = % 0.16g\n",
 586                i, OVI);*/
 587        j++;
 588    
 589        OIII = input[j];
 590        /*fprintf(stderr, "OIII[%d] = % 0.16g\n",
 591                i, OIII);*/
 592        j++;
 593    
 594        OIV = input[j];
 595        /*fprintf(stderr, "OIV[%d] = % 0.16g\n",
 596                i, OIV);*/
 597        j++;
 598    
 599        OV = input[j];
 600        /*fprintf(stderr, "OV[%d] = % 0.16g\n",
 601                i, OV);*/
 602        j++;
 603    
 604        OVIII = input[j];
 605        /*fprintf(stderr, "OVIII[%d] = % 0.16g\n",
 606                i, OVIII);*/
 607        j++;
 608    
 609        OIX = input[j];
 610        /*fprintf(stderr, "OIX[%d] = % 0.16g\n",
 611                i, OIX);*/
 612        j++;
 613
 614        // density = 16*OI + 16*OII + 16*OIII + 16*OIV + 16*OIX + 16*OV + 16*OVI + 16*OVII + 16*OVIII;
 615        /*density = (OI + OII + OIII + OIV + OIX + OV + OVI + OVII + OVIII) * 1e5;
 616        data->Ts[i] = density*ge*mh/(kb*((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))*1e5));// + de/(gamma - 1.0)))*1e5);
 617        data->Ts[i] *= 0.6;*/
 618        data->Ts[i] = ((gamma - 1.0) * 0.6 * mh * ge) / kb; // THIS IS A HACK, BUT IT WILL WORK FOR NOW
 619        if (data->Ts[i] < data->bounds[0]) {
 620            data->Ts[i] = data->bounds[0];
 621        } else if (data->Ts[i] > data->bounds[1]) {
 622            data->Ts[i] = data->bounds[1];
 623        }
 624        data->logTs[i] = log(data->Ts[i]);
 625        /* data->dTs_ge[i] = density*mh/(kb*((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))*1e5));// + de/(gamma - 1.0)))*1e5);
 626         data->dTs_ge[i] *= 0.6;*/
 627        data->dTs_ge[i] = ((gamma - 1.0) * 0.6 * mh) / kb; // THIS IS A HACK, BUT IT WILL WORK FOR NOW
 628        
 629        
 630    
 631        /*fprintf(stderr, "T[%d] = % 0.16g, density = % 0.16g, ge = % 0.16g \n",
 632          i, data->Ts[i], density, ge);*/
 633    
 634    }
 635         
 636}
 637 
 638
 639
 640/*
 641   This setup may be different than the user may anticipate, as a result
 642   of the lockstep timestep we use for a pencil beam through the grid.
 643   As such, it accepts the number of things to interpolate and makes
 644   assumptions about the sizes of the rates.
 645*/
 646
 647/* This also requires no templating other than for the solver name...*/
 648void oxygen_interpolate_rates(oxygen_data *data,
 649                    int nstrip)
 650{
 651    int i, bin_id;
 652    double lb, t1, t2;
 653    lb = log(data->bounds[0]);
 654    /*fprintf(stderr, "lb = % 0.16g, ub = % 0.16g\n", lb, ub);*/
 655    for (i = 0; i < nstrip; i++) {
 656        data->bin_id[i] = bin_id = (int) (data->idbin * (data->logTs[i] - lb));
 657        if (data->bin_id[i] <= 0) {
 658            data->bin_id[i] = 1;
 659        } else if (data->bin_id[i] >= data->nbins) {
 660            data->bin_id[i] = data->nbins - 1;
 661        }
 662        t1 = (lb + (bin_id - 1) * data->dbin);
 663        t2 = (lb + (bin_id    ) * data->dbin);
 664        data->Tdef[i] = (data->logTs[i] - t1)/(t2 - t1);
 665        data->dT[i] = (t2 - t1);
 666        /*fprintf(stderr, "INTERP: %d, bin_id = %d, dT = % 0.16g, T = % 0.16g, logT = % 0.16g\n",
 667                i, data->bin_id[i], data->dT[i], data->Ts[i],
 668                data->logTs[i]);*/
 669    }
 670    for (i = 0; i < nstrip; i++) {
 671        bin_id = data->bin_id[i];
 672        data->rs_OI_i[i] = data->r_OI_i[bin_id] +
 673            data->Tdef[i] * (data->r_OI_i[bin_id+1] - data->r_OI_i[bin_id]);
 674        data->drs_OI_i[i] = (data->r_OI_i[bin_id+1] - data->r_OI_i[bin_id]);
 675        data->drs_OI_i[i] /= data->dT[i];
 676	data->drs_OI_i[i] /= data->Ts[i];
 677    }
 678    
 679    for (i = 0; i < nstrip; i++) {
 680        bin_id = data->bin_id[i];
 681        data->rs_OII_i[i] = data->r_OII_i[bin_id] +
 682            data->Tdef[i] * (data->r_OII_i[bin_id+1] - data->r_OII_i[bin_id]);
 683        data->drs_OII_i[i] = (data->r_OII_i[bin_id+1] - data->r_OII_i[bin_id]);
 684        data->drs_OII_i[i] /= data->dT[i];
 685	data->drs_OII_i[i] /= data->Ts[i];
 686    }
 687    
 688    for (i = 0; i < nstrip; i++) {
 689        bin_id = data->bin_id[i];
 690        data->rs_OII_r[i] = data->r_OII_r[bin_id] +
 691            data->Tdef[i] * (data->r_OII_r[bin_id+1] - data->r_OII_r[bin_id]);
 692        data->drs_OII_r[i] = (data->r_OII_r[bin_id+1] - data->r_OII_r[bin_id]);
 693        data->drs_OII_r[i] /= data->dT[i];
 694	data->drs_OII_r[i] /= data->Ts[i];
 695    }
 696    
 697    for (i = 0; i < nstrip; i++) {
 698        bin_id = data->bin_id[i];
 699        data->rs_OIII_i[i] = data->r_OIII_i[bin_id] +
 700            data->Tdef[i] * (data->r_OIII_i[bin_id+1] - data->r_OIII_i[bin_id]);
 701        data->drs_OIII_i[i] = (data->r_OIII_i[bin_id+1] - data->r_OIII_i[bin_id]);
 702        data->drs_OIII_i[i] /= data->dT[i];
 703	data->drs_OIII_i[i] /= data->Ts[i];
 704    }
 705    
 706    for (i = 0; i < nstrip; i++) {
 707        bin_id = data->bin_id[i];
 708        data->rs_OIII_r[i] = data->r_OIII_r[bin_id] +
 709            data->Tdef[i] * (data->r_OIII_r[bin_id+1] - data->r_OIII_r[bin_id]);
 710        data->drs_OIII_r[i] = (data->r_OIII_r[bin_id+1] - data->r_OIII_r[bin_id]);
 711        data->drs_OIII_r[i] /= data->dT[i];
 712	data->drs_OIII_r[i] /= data->Ts[i];
 713    }
 714    
 715    for (i = 0; i < nstrip; i++) {
 716        bin_id = data->bin_id[i];
 717        data->rs_OIV_i[i] = data->r_OIV_i[bin_id] +
 718            data->Tdef[i] * (data->r_OIV_i[bin_id+1] - data->r_OIV_i[bin_id]);
 719        data->drs_OIV_i[i] = (data->r_OIV_i[bin_id+1] - data->r_OIV_i[bin_id]);
 720        data->drs_OIV_i[i] /= data->dT[i];
 721	data->drs_OIV_i[i] /= data->Ts[i];
 722    }
 723    
 724    for (i = 0; i < nstrip; i++) {
 725        bin_id = data->bin_id[i];
 726        data->rs_OIV_r[i] = data->r_OIV_r[bin_id] +
 727            data->Tdef[i] * (data->r_OIV_r[bin_id+1] - data->r_OIV_r[bin_id]);
 728        data->drs_OIV_r[i] = (data->r_OIV_r[bin_id+1] - data->r_OIV_r[bin_id]);
 729        data->drs_OIV_r[i] /= data->dT[i];
 730	data->drs_OIV_r[i] /= data->Ts[i];
 731    }
 732    
 733    for (i = 0; i < nstrip; i++) {
 734        bin_id = data->bin_id[i];
 735        data->rs_OIX_r[i] = data->r_OIX_r[bin_id] +
 736            data->Tdef[i] * (data->r_OIX_r[bin_id+1] - data->r_OIX_r[bin_id]);
 737        data->drs_OIX_r[i] = (data->r_OIX_r[bin_id+1] - data->r_OIX_r[bin_id]);
 738        data->drs_OIX_r[i] /= data->dT[i];
 739	data->drs_OIX_r[i] /= data->Ts[i];
 740    }
 741    
 742    for (i = 0; i < nstrip; i++) {
 743        bin_id = data->bin_id[i];
 744        data->rs_OV_i[i] = data->r_OV_i[bin_id] +
 745            data->Tdef[i] * (data->r_OV_i[bin_id+1] - data->r_OV_i[bin_id]);
 746        data->drs_OV_i[i] = (data->r_OV_i[bin_id+1] - data->r_OV_i[bin_id]);
 747        data->drs_OV_i[i] /= data->dT[i];
 748	data->drs_OV_i[i] /= data->Ts[i];
 749    }
 750    
 751    for (i = 0; i < nstrip; i++) {
 752        bin_id = data->bin_id[i];
 753        data->rs_OV_r[i] = data->r_OV_r[bin_id] +
 754            data->Tdef[i] * (data->r_OV_r[bin_id+1] - data->r_OV_r[bin_id]);
 755        data->drs_OV_r[i] = (data->r_OV_r[bin_id+1] - data->r_OV_r[bin_id]);
 756        data->drs_OV_r[i] /= data->dT[i];
 757	data->drs_OV_r[i] /= data->Ts[i];
 758    }
 759    
 760    for (i = 0; i < nstrip; i++) {
 761        bin_id = data->bin_id[i];
 762        data->rs_OVI_i[i] = data->r_OVI_i[bin_id] +
 763            data->Tdef[i] * (data->r_OVI_i[bin_id+1] - data->r_OVI_i[bin_id]);
 764        data->drs_OVI_i[i] = (data->r_OVI_i[bin_id+1] - data->r_OVI_i[bin_id]);
 765        data->drs_OVI_i[i] /= data->dT[i];
 766	data->drs_OVI_i[i] /= data->Ts[i];
 767    }
 768    
 769    for (i = 0; i < nstrip; i++) {
 770        bin_id = data->bin_id[i];
 771        data->rs_OVI_r[i] = data->r_OVI_r[bin_id] +
 772            data->Tdef[i] * (data->r_OVI_r[bin_id+1] - data->r_OVI_r[bin_id]);
 773        data->drs_OVI_r[i] = (data->r_OVI_r[bin_id+1] - data->r_OVI_r[bin_id]);
 774        data->drs_OVI_r[i] /= data->dT[i];
 775	data->drs_OVI_r[i] /= data->Ts[i];
 776    }
 777    
 778    for (i = 0; i < nstrip; i++) {
 779        bin_id = data->bin_id[i];
 780        data->rs_OVII_i[i] = data->r_OVII_i[bin_id] +
 781            data->Tdef[i] * (data->r_OVII_i[bin_id+1] - data->r_OVII_i[bin_id]);
 782        data->drs_OVII_i[i] = (data->r_OVII_i[bin_id+1] - data->r_OVII_i[bin_id]);
 783        data->drs_OVII_i[i] /= data->dT[i];
 784	data->drs_OVII_i[i] /= data->Ts[i];
 785    }
 786    
 787    for (i = 0; i < nstrip; i++) {
 788        bin_id = data->bin_id[i];
 789        data->rs_OVII_r[i] = data->r_OVII_r[bin_id] +
 790            data->Tdef[i] * (data->r_OVII_r[bin_id+1] - data->r_OVII_r[bin_id]);
 791        data->drs_OVII_r[i] = (data->r_OVII_r[bin_id+1] - data->r_OVII_r[bin_id]);
 792        data->drs_OVII_r[i] /= data->dT[i];
 793	data->drs_OVII_r[i] /= data->Ts[i];
 794    }
 795    
 796    for (i = 0; i < nstrip; i++) {
 797        bin_id = data->bin_id[i];
 798        data->rs_OVIII_i[i] = data->r_OVIII_i[bin_id] +
 799            data->Tdef[i] * (data->r_OVIII_i[bin_id+1] - data->r_OVIII_i[bin_id]);
 800        data->drs_OVIII_i[i] = (data->r_OVIII_i[bin_id+1] - data->r_OVIII_i[bin_id]);
 801        data->drs_OVIII_i[i] /= data->dT[i];
 802	data->drs_OVIII_i[i] /= data->Ts[i];
 803    }
 804    
 805    for (i = 0; i < nstrip; i++) {
 806        bin_id = data->bin_id[i];
 807        data->rs_OVIII_r[i] = data->r_OVIII_r[bin_id] +
 808            data->Tdef[i] * (data->r_OVIII_r[bin_id+1] - data->r_OVIII_r[bin_id]);
 809        data->drs_OVIII_r[i] = (data->r_OVIII_r[bin_id+1] - data->r_OVIII_r[bin_id]);
 810        data->drs_OVIII_r[i] /= data->dT[i];
 811	data->drs_OVIII_r[i] /= data->Ts[i];
 812    }
 813    
 814
 815}
 816 
 817
 818
 819
 820int calculate_rhs_oxygen(double *input, double *rhs, int nstrip,
 821                  int nchem, void *sdata)
 822{
 823    /* We iterate over all of the rates */
 824    /* Calculate temperature first */
 825    oxygen_data *data = (oxygen_data*)sdata;
 826    int i, j;
 827    oxygen_calculate_temperature(data, input, nstrip, nchem);
 828
 829    oxygen_interpolate_rates(data, nstrip);
 830
 831    /* Now we set up some temporaries */
 832    //double *T = data->Ts;
 833    double *OI_i = data->rs_OI_i;
 834    double *OII_i = data->rs_OII_i;
 835    double *OII_r = data->rs_OII_r;
 836    double *OIII_i = data->rs_OIII_i;
 837    double *OIII_r = data->rs_OIII_r;
 838    double *OIV_i = data->rs_OIV_i;
 839    double *OIV_r = data->rs_OIV_r;
 840    double *OIX_r = data->rs_OIX_r;
 841    double *OV_i = data->rs_OV_i;
 842    double *OV_r = data->rs_OV_r;
 843    double *OVI_i = data->rs_OVI_i;
 844    double *OVI_r = data->rs_OVI_r;
 845    double *OVII_i = data->rs_OVII_i;
 846    double *OVII_r = data->rs_OVII_r;
 847    double *OVIII_i = data->rs_OVIII_i;
 848    double *OVIII_r = data->rs_OVIII_r;
 849    double ge;
 850    double de;
 851    double OI;
 852    double OVII;
 853    double OII;
 854    double OVI;
 855    double OIII;
 856    double OIV;
 857    double OV;
 858    double OVIII;
 859    double OIX;
 860
 861    double mh = 1.67e-24;
 862    double total, total_e, total_de, mdensity;
 863    double floor_value = 1e-25; //this should probably be globally stored
 864    int is_floored[nchem];
 865    for (i = 0; i<nstrip; i++) {
 866      //fprintf(stderr, "T = %0.5g\n", T[i]);
 867        
 868        for (int s = 0; s < nchem; s++)
 869            is_floored[s] = 0;
 870        
 871        j = i * nchem;
 872        total = total_e = total_de = mdensity = 0.0;
 873        ge = input[j];
 874        if (ge < 0.0) {
 875          //fprintf(stderr, "RNegative[%d][ge] = % 0.16g [%d]\n",
 876          //  i, ge, j);
 877            return 1;
 878          ge = 1e-20;
 879        }
 880        if (ge <= floor_value) {
 881            is_floored[j % nchem] = 1;
 882        }
 883
 884        
 885        j++;
 886    
 887        de = input[j];
 888        if (de < 0.0) {
 889          //fprintf(stderr, "RNegative[%d][de] = % 0.16g [%d]\n",
 890          //  i, de, j);
 891            return 1;
 892          de = 1e-20;
 893        }
 894	if (de <= floor_value) {
 895          is_floored[j % nchem] = 1;
 896        }
 897        
 898
 899        j++;
 900    
 901        OI = input[j];
 902        if (OI < 0.0) {
 903          //fprintf(stderr, "RNegative[%d][OI] = % 0.16g [%d]\n",
 904          //  i, OI, j);
 905            return 1;
 906          OI = 1e-20;
 907        }
 908	if (OI <= floor_value) {
 909          is_floored[j % nchem] = 1;
 910        }        
 911        
 912          total+=OI * 16;
 913        
 914        
 915        j++;
 916    
 917        OVII = input[j];
 918        if (OVII < 0.0) {
 919          //fprintf(stderr, "RNegative[%d][OVII] = % 0.16g [%d]\n",
 920          //  i, OVII, j);
 921            return 1;
 922          OVII = 1e-20;
 923        }
 924	if (OVII <= floor_value) {
 925          is_floored[j % nchem] = 1;
 926        }        
 927        
 928          total+=OVII * 16;
 929        
 930        
 931        j++;
 932    
 933        OII = input[j];
 934        if (OII < 0.0) {
 935          //fprintf(stderr, "RNegative[%d][OII] = % 0.16g [%d]\n",
 936          //  i, OII, j);
 937            return 1;
 938          OII = 1e-20;
 939        }
 940	if (OII <= floor_value) {
 941          is_floored[j % nchem] = 1;
 942        }
 943        
 944          total+=OII * 16;
 945        
 946        
 947        j++;
 948    
 949        OVI = input[j];
 950        if (OVI < 0.0) {
 951          //fprintf(stderr, "RNegative[%d][OVI] = % 0.16g [%d]\n",
 952          //  i, OVI, j);
 953            return 1;
 954          OVI = 1e-20;
 955        }
 956	if (OVI <= floor_value) {
 957          is_floored[j % nchem] = 1;
 958        }
 959        
 960          total+=OVI * 16;
 961        
 962        
 963        j++;
 964    
 965        OIII = input[j];
 966        if (OIII < 0.0) {
 967          //fprintf(stderr, "RNegative[%d][OIII] = % 0.16g [%d]\n",
 968          //  i, OIII, j);
 969            return 1;
 970          OIII = 1e-20;
 971        }
 972	if (OIII <= floor_value) {
 973          is_floored[j % nchem] = 1;
 974        }
 975        
 976          total+=OIII * 16;
 977        
 978        
 979        j++;
 980    
 981        OIV = input[j];
 982        if (OIV < 0.0) {
 983          //fprintf(stderr, "RNegative[%d][OIV] = % 0.16g [%d]\n",
 984          //  i, OIV, j);
 985            return 1;
 986          OIV = 1e-20;
 987        }
 988	if (OIV <= floor_value) {
 989          is_floored[j % nchem] = 1;
 990        }
 991        
 992          total+=OIV * 16;
 993        
 994        
 995        j++;
 996    
 997        OV = input[j];
 998        if (OV < 0.0) {
 999          //fprintf(stderr, "RNegative[%d][OV] = % 0.16g [%d]\n",
1000          //  i, OV, j);
1001            return 1;
1002          OV = 1e-20;
1003        }
1004	if (OV <= floor_value) {
1005          is_floored[j % nchem] = 1;
1006        }
1007        
1008        
1009          total+=OV * 16;
1010        
1011        
1012        j++;
1013    
1014        OVIII = input[j];
1015        if (OVIII < 0.0) {
1016          //fprintf(stderr, "RNegative[%d][OVIII] = % 0.16g [%d]\n",
1017          //  i, OVIII, j);
1018            return 1;
1019          OVIII = 1e-20;
1020        }
1021	if (OVIII <= floor_value) {
1022          is_floored[j % nchem] = 1;
1023        }
1024        
1025          total+=OVIII * 16;
1026        
1027        
1028        j++;
1029    
1030        OIX = input[j];
1031        if (OIX < 0.0) {
1032          //fprintf(stderr, "RNegative[%d][OIX] = % 0.16g [%d]\n",
1033          //  i, OIX, j);
1034            return 1;
1035          OIX = 1e-20;
1036        }
1037	if (OIX <= floor_value) {
1038          is_floored[j % nchem] = 1;
1039        }
1040        
1041          total+=OIX * 16;
1042        
1043        
1044        j++;
1045    
1046        mdensity = total * mh;
1047        total = 0.0;
1048        j = i * nchem;
1049        // 
1050        // Species: ge
1051        // 
1052        rhs[j] = 0.0;
1053        
1054	//rhs[j] /= mdensity;
1055        
1056        
1057        j++;
1058    
1059        // 
1060        // Species: de
1061        // 
1062        rhs[j] = OIII_i[i]*OIII*de - OIII_r[i]*OIII*de + OII_i[i]*OII*de - OII_r[i]*OII*de + OIV_i[i]*OIV*de - OIV_r[i]*OIV*de - OIX_r[i]*OIX*de + OI_i[i]*OI*de + OVIII_i[i]*OVIII*de - OVIII_r[i]*OVIII*de + OVII_i[i]*OVII*de - OVII_r[i]*OVII*de + OVI_i[i]*OVI*de - OVI_r[i]*OVI*de + OV_i[i]*OV*de - OV_r[i]*OV*de;
1063        
1064            total_de += -rhs[j];
1065        
1066        j++;
1067    
1068        // 
1069        // Species: OI
1070        // 
1071        rhs[j] = OII_r[i]*OII*de - OI_i[i]*OI*de;
1072
1073            /* Already in number density, not mass density */
1074            total += rhs[j] * 16;
1075            total_e += OI * 0;
1076        
1077        
1078            total_de += rhs[j] * 0;
1079        
1080        j++;
1081    
1082        // 
1083        // Species: OVII
1084        // 
1085        rhs[j] = OVIII_r[i]*OVIII*de - OVII_i[i]*OVII*de - OVII_r[i]*OVII*de + OVI_i[i]*OVI*de;
1086
1087            /* Already in number density, not mass density */
1088            total += rhs[j] * 16;
1089            total_e += OVII * 6;
1090        
1091        
1092            total_de += rhs[j] * 6;
1093        
1094        j++;
1095    
1096        // 
1097        // Species: OII
1098        // 
1099        rhs[j] = OIII_r[i]*OIII*de - OII_i[i]*OII*de - OII_r[i]*OII*de + OI_i[i]*OI*de;
1100
1101            /* Already in number density, not mass density */
1102            total += rhs[j] * 16;
1103            total_e += OII * 1;
1104        
1105        
1106            total_de += rhs[j] * 1;
1107        
1108        j++;
1109    
1110        // 
1111        // Species: OVI
1112        // 
1113        rhs[j] = OVII_r[i]*OVII*de - OVI_i[i]*OVI*de - OVI_r[i]*OVI*de + OV_i[i]*OV*de;
1114
1115            /* Already in number density, not mass density */
1116            total += rhs[j] * 16;
1117            total_e += OVI * 5;
1118        
1119        
1120            total_de += rhs[j] * 5;
1121        
1122        j++;
1123    
1124        // 
1125        // Species: OIII
1126        // 
1127        rhs[j] = -OIII_i[i]*OIII*de - OIII_r[i]*OIII*de + OII_i[i]*OII*de + OIV_r[i]*OIV*de;
1128
1129            /* Already in number density, not mass density */
1130            total += rhs[j] * 16;
1131            total_e += OIII * 2;
1132        
1133        
1134            total_de += rhs[j] * 2;
1135        
1136        j++;
1137    
1138        // 
1139        // Species: OIV
1140        // 
1141        rhs[j] = OIII_i[i]*OIII*de - OIV_i[i]*OIV*de - OIV_r[i]*OIV*de + OV_r[i]*OV*de;
1142
1143            /* Already in number density, not mass density */
1144            total += rhs[j] * 16;
1145            total_e += OIV * 3;
1146        
1147        
1148            total_de += rhs[j] * 3;
1149        
1150        j++;
1151    
1152        // 
1153        // Species: OV
1154        // 
1155        rhs[j] = OIV_i[i]*OIV*de + OVI_r[i]*OVI*de - OV_i[i]*OV*de - OV_r[i]*OV*de;
1156
1157            /* Already in number density, not mass density */
1158            total += rhs[j] * 16;
1159            total_e += OV * 4;
1160        
1161        
1162            total_de += rhs[j] * 4;
1163        
1164        j++;
1165    
1166        // 
1167        // Species: OVIII
1168        // 
1169        rhs[j] = OIX_r[i]*OIX*de - OVIII_i[i]*OVIII*de - OVIII_r[i]*OVIII*de + OVII_i[i]*OVII*de;
1170
1171            /* Already in number density, not mass density */
1172            total += rhs[j] * 16;
1173            total_e += OVIII * 7;
1174        
1175        
1176            total_de += rhs[j] * 7;
1177        
1178        j++;
1179    
1180        // 
1181        // Species: OIX
1182        // 
1183        rhs[j] = -OIX_r[i]*OIX*de + OVIII_i[i]*OVIII*de;
1184
1185            /* Already in number density, not mass density */
1186            total += rhs[j] * 16;
1187            total_e += OIX * 8;
1188        
1189        
1190            total_de += rhs[j] * 8;
1191        
1192        j++;
1193    
1194    }  
1195    return 0;
1196}
1197
1198
1199
1200
1201int calculate_jacobian_oxygen(double *input, double *Joutput,
1202        int nstrip, int nchem, void *sdata)
1203{
1204    /* We iterate over all of the rates */
1205    /* Calculate temperature first */
1206    oxygen_data *data = (oxygen_data*)sdata;
1207
1208    int i, j;
1209    oxygen_calculate_temperature(data, input, nstrip, nchem);
1210
1211    oxygen_interpolate_rates(data, nstrip);
1212
1213    /* Now we set up some temporaries */
1214    double *Tge = data->dTs_ge;
1215    double *OI_i = data->rs_OI_i;
1216    double *rOI_i = data->drs_OI_i;
1217    double *OII_i = data->rs_OII_i;
1218    double *rOII_i = data->drs_OII_i;
1219    double *OII_r = data->rs_OII_r;
1220    double *rOII_r = data->drs_OII_r;
1221    double *OIII_i = data->rs_OIII_i;
1222    double *rOIII_i = data->drs_OIII_i;
1223    double *OIII_r = data->rs_OIII_r;
1224    double *rOIII_r = data->drs_OIII_r;
1225    double *OIV_i = data->rs_OIV_i;
1226    double *rOIV_i = data->drs_OIV_i;
1227    double *OIV_r = data->rs_OIV_r;
1228    double *rOIV_r = data->drs_OIV_r;
1229    double *OIX_r = data->rs_OIX_r;
1230    double *rOIX_r = data->drs_OIX_r;
1231    double *OV_i = data->rs_OV_i;
1232    double *rOV_i = data->drs_OV_i;
1233    double *OV_r = data->rs_OV_r;
1234    double *rOV_r = data->drs_OV_r;
1235    double *OVI_i = data->rs_OVI_i;
1236    double *rOVI_i = data->drs_OVI_i;
1237    double *OVI_r = data->rs_OVI_r;
1238    double *rOVI_r = data->drs_OVI_r;
1239    double *OVII_i = data->rs_OVII_i;
1240    double *rOVII_i = data->drs_OVII_i;
1241    double *OVII_r = data->rs_OVII_r;
1242    double *rOVII_r = data->drs_OVII_r;
1243    double *OVIII_i = data->rs_OVIII_i;
1244    double *rOVIII_i = data->drs_OVIII_i;
1245    double *OVIII_r = data->rs_OVIII_r;
1246    double *rOVIII_r = data->drs_OVIII_r;
1247    double ge;
1248    double de;
1249    double OI;
1250    double OVII;
1251    double OII;
1252    double OVI;
1253    double OIII;
1254    double OIV;
1255    double OV;
1256    double OVIII;
1257    double OIX;
1258
1259    double mh = 1.67e-24;
1260    double total, mdensity;
1261    for (i = 0; i<nstrip; i++) {
1262        j = i * nchem;
1263	total = mdensity = 0.0;
1264	    ge = input[j];
1265        if (ge < 0.0) {
1266          fprintf(stderr, "JNegative[%d][ge] = % 0.16g [%d]\n",
1267            i, ge, j);
1268          /*ge = 0.0;*/
1269          ge = 1e-20;
1270          return 1;
1271        }
1272	
1273        j++;
1274    
1275	    de = input[j];
1276        if (de < 0.0) {
1277          fprintf(stderr, "JNegative[%d][de] = % 0.16g [%d]\n",
1278            i, de, j);
1279          /*de = 0.0;*/
1280          de = 1e-20;
1281          return 1;
1282        }
1283	
1284        
1285	
1286        j++;
1287    
1288	    OI = input[j];
1289        if (OI < 0.0) {
1290          fprintf(stderr, "JNegative[%d][OI] = % 0.16g [%d]\n",
1291            i, OI, j);
1292          /*OI = 0.0;*/
1293          OI = 1e-20;
1294          return 1;
1295        }
1296	
1297        
1298          total+=OI * 16;
1299        
1300	
1301        j++;
1302    
1303	    OVII = input[j];
1304        if (OVII < 0.0) {
1305          fprintf(stderr, "JNegative[%d][OVII] = % 0.16g [%d]\n",
1306            i, OVII, j);
1307          /*OVII = 0.0;*/
1308          OVII = 1e-20;
1309          return 1;
1310        }
1311	
1312        
1313          total+=OVII * 16;
1314        
1315	
1316        j++;
1317    
1318	    OII = input[j];
1319        if (OII < 0.0) {
1320          fprintf(stderr, "JNegative[%d][OII] = % 0.16g [%d]\n",
1321            i, OII, j);
1322          /*OII = 0.0;*/
1323          OII = 1e-20;
1324          return 1;
1325        }
1326	
1327        
1328          total+=OII * 16;
1329        
1330	
1331        j++;
1332    
1333	    OVI = input[j];
1334        if (OVI < 0.0) {
1335          fprintf(stderr, "JNegative[%d][OVI] = % 0.16g [%d]\n",
1336            i, OVI, j);
1337          /*OVI = 0.0;*/
1338          OVI = 1e-20;
1339          return 1;
1340        }
1341	
1342        
1343          total+=OVI * 16;
1344        
1345	
1346        j++;
1347    
1348	    OIII = input[j];
1349        if (OIII < 0.0) {
1350          fprintf(stderr, "JNegative[%d][OIII] = % 0.16g [%d]\n",
1351            i, OIII, j);
1352          /*OIII = 0.0;*/
1353          OIII = 1e-20;
1354          return 1;
1355        }
1356	
1357        
1358          total+=OIII * 16;
1359        
1360	
1361        j++;
1362    
1363	    OIV = input[j];
1364        if (OIV < 0.0) {
1365          fprintf(stderr, "JNegative[%d][OIV] = % 0.16g [%d]\n",
1366            i, OIV, j);
1367          /*OIV = 0.0;*/
1368          OIV = 1e-20;
1369          return 1;
1370        }
1371	
1372        
1373          total+=OIV * 16;
1374        
1375	
1376        j++;
1377    
1378	    OV = input[j];
1379        if (OV < 0.0) {
1380          fprintf(stderr, "JNegative[%d][OV] = % 0.16g [%d]\n",
1381            i, OV, j);
1382          /*OV = 0.0;*/
1383          OV = 1e-20;
1384          return 1;
1385        }
1386	
1387        
1388          total+=OV * 16;
1389        
1390	
1391        j++;
1392    
1393	    OVIII = input[j];
1394        if (OVIII < 0.0) {
1395          fprintf(stderr, "JNegative[%d][OVIII] = % 0.16g [%d]\n",
1396            i, OVIII, j);
1397          /*OVIII = 0.0;*/
1398          OVIII = 1e-20;
1399          return 1;
1400        }
1401	
1402        
1403          total+=OVIII * 16;
1404        
1405	
1406        j++;
1407    
1408	    OIX = input[j];
1409        if (OIX < 0.0) {
1410          fprintf(stderr, "JNegative[%d][OIX] = % 0.16g [%d]\n",
1411            i, OIX, j);
1412          /*OIX = 0.0;*/
1413          OIX = 1e-20;
1414          return 1;
1415        }
1416	
1417        
1418          total+=OIX * 16;
1419        
1420	
1421        j++;
1422    
1423        mdensity = total * mh;
1424        
1425        j = i * nchem * nchem;
1426        // 
1427        // Species: ge
1428        //
1429            // ge by ge
1430            Joutput[j] = 0.0;
1431	    
1432	    //Joutput[j] /= mdensity;
1433	    
1434	    
1435            Joutput[j] *= Tge[i];
1436            
1437            j++;
1438            // de by ge
1439            Joutput[j] = OI*de*rOI_i[i] + OII*de*rOII_i[i] - OII*de*rOII_r[i] + OIII*de*rOIII_i[i] - OIII*de*rOIII_r[i] + OIV*de*rOIV_i[i] - OIV*de*rOIV_r[i] - OIX*de*rOIX_r[i] + OV*de*rOV_i[i] - OV*de*rOV_r[i] + OVI*de*rOVI_i[i] - OVI*de*rOVI_r[i] + OVII*de*rOVII_i[i] - OVII*de*rOVII_r[i] + OVIII*de*rOVIII_i[i] - OVIII*de*rOVIII_r[i];
1440	    
1441	    
1442            Joutput[j] *= Tge[i];
1443            
1444            j++;
1445            // OI by ge
1446            Joutput[j] = -OI*de*rOI_i[i] + OII*de*rOII_r[i];
1447	    
1448	    
1449            Joutput[j] *= Tge[i];
1450            
1451            j++;
1452            // OVII by ge
1453            Joutput[j] = OVI*de*rOVI_i[i] - OVII*de*rOVII_i[i] - OVII*de*rOVII_r[i] + OVIII*de*rOVIII_r[i];
1454	    
1455	    
1456            Joutput[j] *= Tge[i];
1457            
1458            j++;
1459            // OII by ge
1460            Joutput[j] = OI*de*rOI_i[i] - OII*de*rOII_i[i] - OII*de*rOII_r[i] + OIII*de*rOIII_r[i];
1461	    
1462	    
1463            Joutput[j] *= Tge[i];
1464            
1465            j++;
1466            // OVI by ge
1467            Joutput[j] = OV*de*rOV_i[i] - OVI*de*rOVI_i[i] - OVI*de*rOVI_r[i] + OVII*de*rOVII_r[i];
1468	    
1469	    
1470            Joutput[j] *= Tge[i];
1471            
1472            j++;
1473            // OIII by ge
1474            Joutput[j] = OII*de*rOII_i[i] - OIII*de*rOIII_i[i] - OIII*de*rOIII_r[i] + OIV*de*rOIV_r[i];
1475	    
1476	    
1477            Joutput[j] *= Tge[i];
1478            
1479            j++;
1480            // OIV by ge
1481            Joutput[j] = OIII*de*rOIII_i[i] - OIV*de*rOIV_i[i] - OIV*de*rOIV_r[i] + OV*de*rOV_r[i];
1482	    
1483	    
1484            Joutput[j] *= Tge[i];
1485            
1486            j++;
1487            // OV by ge
1488            Joutput[j] = OIV*de*rOIV_i[i] - OV*de*rOV_i[i] - OV*de*rOV_r[i] + OVI*de*rOVI_r[i];
1489	    
1490	    
1491            Joutput[j] *= Tge[i];
1492            
1493            j++;
1494            // OVIII by ge
1495            Joutput[j] = OIX*de*rOIX_r[i] + OVII*de*rOVII_i[i] - OVIII*de*rOVIII_i[i] - OVIII*de*rOVIII_r[i];
1496	    
1497	    
1498            Joutput[j] *= Tge[i];
1499            
1500            j++;
1501            // OIX by ge
1502            Joutput[j] = -OIX*de*rOIX_r[i] + OVIII*de*rOVIII_i[i];
1503	    
1504	    
1505            Joutput[j] *= Tge[i];
1506            
1507            j++;
1508    
1509        // 
1510        // Species: de
1511        //
1512            // ge by de
1513            Joutput[j] = 0.0;
1514	    
1515	    //Joutput[j] /= mdensity;
1516	    
1517	    
1518            j++;
1519            // de by de
1520            Joutput[j] = 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;
1521	    
1522	    
1523            j++;
1524            // OI by de
1525            Joutput[j] = OII_r[i]*OII - OI_i[i]*OI;
1526	    
1527	    
1528            j++;
1529            // OVII by de
1530            Joutput[j] = OVIII_r[i]*OVIII - OVII_i[i]*OVII - OVII_r[i]*OVII + OVI_i[i]*OVI;
1531	    
1532	    
1533            j++;
1534            // OII by de
1535            Joutput[j] = OIII_r[i]*OIII - OII_i[i]*OII - OII_r[i]*OII + OI_i[i]*OI;
1536	    
1537	    
1538            j++;
1539            // OVI by de
1540            Joutput[j] = OVII_r[i]*OVII - OVI_i[i]*OVI - OVI_r[i]*OVI + OV_i[i]*OV;
1541	    
1542	    
1543            j++;
1544            // OIII by de
1545            Joutput[j] = -OIII_i[i]*OIII - OIII_r[i]*OIII + OII_i[i]*OII + OIV_r[i]*OIV;
1546	    
1547	    
1548            j++;
1549            // OIV by de
1550            Joutput[j] = OIII_i[i]*OIII - OIV_i[i]*OIV - OIV_r[i]*OIV + OV_r[i]*OV;
1551	    
1552	    
1553            j++;
1554            // OV by de
1555            Joutput[j] = OIV_i[i]*OIV + OVI_r[i]*OVI - OV_i[i]*OV - OV_r[i]*OV;
1556	    
1557	    
1558            j++;
1559            // OVIII by de
1560            Joutput[j] = OIX_r[i]*OIX - OVIII_i[i]*OVIII - OVIII_r[i]*OVIII + OVII_i[i]*OVII;
1561	    
1562	    
1563            j++;
1564            // OIX by de
1565            Joutput[j] = -OIX_r[i]*OIX + OVIII_i[i]*OVIII;
1566	    
1567	    
1568            j++;
1569    
1570        // 
1571        // Species: OI
1572        //
1573            // ge by OI
1574            Joutput[j] = 0.0;
1575	    
1576	    //Joutput[j] /= mdensity;
1577	    
1578	    
1579            j++;
1580            // de by OI
1581            Joutput[j] = OI_i[i]*de;
1582	    
1583	    
1584            j++;
1585            // OI by OI
1586            Joutput[j] = -OI_i[i]*de;
1587	    
1588	    
1589            j++;
1590            // OVII by OI
1591            Joutput[j] = 0.0;
1592	    
1593	    
1594            j++;
1595            // OII by OI
1596            Joutput[j] = OI_i[i]*de;
1597	    
1598	    
1599            j++;
1600            // OVI by OI
1601            Joutput[j] = 0.0;
1602	    
1603	    
1604            j++;
1605            // OIII by OI
1606            Joutput[j] = 0.0;
1607	    
1608	    
1609            j++;
1610            // OIV by OI
1611            Joutput[j] = 0.0;
1612	    
1613	    
1614            j++;
1615            // OV by OI
1616            Joutput[j] = 0.0;
1617	    
1618	    
1619            j++;
1620            // OVIII by OI
1621            Joutput[j] = 0.0;
1622	    
1623	    
1624            j++;
1625            // OIX by OI
1626            Joutput[j] = 0.0;
1627	    
1628	    
1629            j++;
1630    
1631        // 
1632        // Species: OVII
1633        //
1634            // ge by OVII
1635            Joutput[j] = 0.0;
1636	    
1637	    //Joutput[j] /= mdensity;
1638	    
1639	    
1640            j++;
1641            // de by OVII
1642            Joutput[j] = OVII_i[i]*de - OVII_r[i]*de;
1643	    
1644	    
1645            j++;
1646            // OI by OVII
1647            Joutput[j] = 0.0;
1648	    
1649	    
1650            j++;
1651            // OVII by OVII
1652            Joutput[j] = -OVII_i[i]*de - OVII_r[i]*de;
1653	    
1654	    
1655            j++;
1656            // OII by OVII
1657            Joutput[j] = 0.0;
1658	    
1659	    
1660            j++;
1661            // OVI by OVII
1662            Joutput[j] = OVII_r[i]*de;
1663	    
1664	    
1665            j++;
1666            // OIII by OVII
1667            Joutput[j] = 0.0;
1668	    
1669	    
1670            j++;
1671            // OIV by OVII
1672            Joutput[j] = 0.0;
1673	    
1674	    
1675            j++;
1676            // OV by OVII
1677            Joutput[j] = 0.0;
1678	    
1679	    
1680            j++;
1681            // OVIII by OVII
1682            Joutput[j] = OVII_i[i]*de;
1683	    
1684	    
1685            j++;
1686            // OIX by OVII
1687            Joutput[j] = 0.0;
1688	    
1689	    
1690            j++;
1691    
1692        // 
1693        // Species: OII
1694        //
1695            // ge by OII
1696            Joutput[j] = 0.0;
1697	    
1698	    //Joutput[j] /= mdensity;
1699	    
1700	    
1701            j++;
1702            // de by OII
1703            Joutput[j] = OII_i[i]*de - OII_r[i]*de;
1704	    
1705	    
1706            j++;
1707            // OI by OII
1708            Joutput[j] = OII_r[i]*de;
1709	    
1710	    
1711            j++;
1712            // OVII by OII
1713            Joutput[j] = 0.0;
1714	    
1715	    
1716            j++;
1717            // OII by OII
1718            Joutput[j] = -OII_i[i]*de - OII_r[i]*de;
1719	    
1720	    
1721            j++;
1722            // OVI by OII
1723            Joutput[j] = 0.0;
1724	    
1725	    
1726            j++;
1727            // OIII by OII
1728            Joutput[j] = OII_i[i]*de;
1729	    
1730	    
1731            j++;
1732            // OIV by OII
1733            Joutput[j] = 0.0;
1734	    
1735	    
1736            j++;
1737            // OV by OII
1738            Joutput[j] = 0.0;
1739	    
1740	    
1741            j++;
1742            // OVIII by OII
1743            Joutput[j] = 0.0;
1744	    
1745	    
1746            j++;
1747            // OIX by OII
1748            Joutput[j] = 0.0;
1749	    
1750	    
1751            j++;
1752    
1753        // 
1754        // Species: OVI
1755        //
1756            // ge by OVI
1757            Joutput[j] = 0.0;
1758	    
1759	    //Joutput[j] /= mdensity;
1760	    
1761	    
1762            j++;
1763            // de by OVI
1764            Joutput[j] = OVI_i[i]*de - OVI_r[i]*de;
1765	    
1766	    
1767            j++;
1768            // OI by OVI
1769            Joutput[j] = 0.0;
1770	    
1771	    
1772            j++;
1773            // OVII by OVI
1774            Joutput[j] = OVI_i[i]*de;
1775	    
1776	    
1777            j++;
1778            // OII by OVI
1779            Joutput[j] = 0.0;
1780	    
1781	    
1782            j++;
1783            // OVI by OVI
1784            Joutput[j] = -OVI_i[i]*de - OVI_r[i]*de;
1785	    
1786	    
1787            j++;
1788            

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