diff --git a/src/MINT/MINT_init_store.c b/src/MINT/MINT_init_store.c index 91ec512d66f41b9b2e0a14da89df66c0bb7d9186..a52495740d8d3acee399250d2c259fec2cb38dac 100644 --- a/src/MINT/MINT_init_store.c +++ b/src/MINT/MINT_init_store.c @@ -14,6 +14,7 @@ void MINT_init_store(struct store_t * store) store->MINT_tables[i] = NULL; store->MINT_data_available[i] = NULL; store->MINT_data_map[i] = NULL; + store->MINT_initial_coordinate[i] = -1.0; } store->MINT_loaded = FALSE; } diff --git a/src/MINT/MINT_initial_XHc.c b/src/MINT/MINT_initial_XHc.c index e83120ee75c2f311f1186baf6c89cee53ee7b36b..938ec498378e1f20df9053ad2e217c23d76aff79 100644 --- a/src/MINT/MINT_initial_XHc.c +++ b/src/MINT/MINT_initial_XHc.c @@ -4,40 +4,49 @@ No_empty_translation_unit_warning; #ifdef MINT #include "MINT.h" - double MINT_initial_XHc(struct stardata_t * const stardata) { /* * Determine the initial central hydrogen abundance * in the MINT models */ - /* - const double MINT_offset = 0.7 - 0.6985; - const double Z = Max(0.0, - 0.760 - 3.0 * stardata->common.metallicity - MINT_offset); - return Z; - */ + if(stardata->store->MINT_initial_coordinate[MINT_TABLE_MS] < -TINY) + { + const size_t nperline = + stardata->store->MINT_tables[MINT_TABLE_MS]->ndata + + stardata->store->MINT_tables[MINT_TABLE_MS]->nparam; + double * const pmin = stardata->store->MINT_tables[MINT_TABLE_MS]->data + MINT_MS_CENTRAL_HYDROGEN; + double * const pmax = pmin + nperline * (1 + stardata->store->MINT_tables[MINT_TABLE_MS]->nlines); - int i = 1; // XHC column - const int nperline = - stardata->store->MINT_tables[MINT_TABLE_MS]->ndata + - stardata->store->MINT_tables[MINT_TABLE_MS]->nparam; + /* + * While the abundance is increasing, keep looping through + * lines of data. XHc will increase, but then when you + * get to the next star (or the end of the table) it will decrease. + * + * The pointer p points to the hydrogen abundance at each line. + * You should increment this by nperline. + * + * Note that pmax is the last value +1 line because p<pmax (not <=). + */ + double * p = pmin; + double prev_XHC0 = *p; + while(p < pmax && *p >= prev_XHC0) + { + prev_XHC0 = *p; + p += nperline; + } - double * p = stardata->store->MINT_tables[MINT_TABLE_MS]->data + i; - double * const pmin = p; - double * const pmax = p + nperline * stardata->store->MINT_tables[MINT_TABLE_MS]->nlines; + /* we've overshot */ + p -= nperline; - const double XHC0 = *p; /* first value */ + /* make sure we haven't wrapped back around the first value */ + p = Max(p, pmin); - /* - * While the abundance is increasing, keep looping through - * lines of data. When it - */ - while(p < pmax && *p >= XHC0) - { - p += nperline; + Dprint("Initial hydrogen abundance %g\n",*p); + stardata->store->MINT_initial_coordinate[MINT_TABLE_MS] = *p; } - p = Max(p - nperline, pmin); - return *p; + + /* hence return the central hydrogen */ + return stardata->store->MINT_initial_coordinate[MINT_TABLE_MS]; } #endif //MINT diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h index 048b1823e566cf8bae66c9d13906b5089ecb308e..7b86de4e46d034946306a38700b5db3f9dabb0c7 100644 --- a/src/binary_c_structures.h +++ b/src/binary_c_structures.h @@ -363,6 +363,7 @@ struct store_t { struct mint_header_t * MINT_headers[MINT_TABLE_NUMBER]; unsigned int * MINT_data_map[MINT_TABLE_NUMBER]; unsigned int MINT_tables_nscalars[MINT_TABLE_NUMBER]; + Abundance MINT_initial_coordinate[MINT_TABLE_NUMBER]; #endif//MINT int * mint_generic_map[NUMBER_OF_STELLAR_TYPES];