From 60ada3b7e33c2cbe4ed6c5f1e393c052a0ca57ae Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Sat, 7 Aug 2021 15:47:06 +0100
Subject: [PATCH] cache the initial XHc

---
 src/MINT/MINT_init_store.c  |  1 +
 src/MINT/MINT_initial_XHc.c | 57 +++++++++++++++++++++----------------
 src/binary_c_structures.h   |  1 +
 3 files changed, 35 insertions(+), 24 deletions(-)

diff --git a/src/MINT/MINT_init_store.c b/src/MINT/MINT_init_store.c
index 91ec512d6..a52495740 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 e83120ee7..938ec4983 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 048b1823e..7b86de4e4 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];
 
-- 
GitLab