From e4a6f314b56050486afadde3ced0552dbc253212 Mon Sep 17 00:00:00 2001
From: Yago Ascasibar <yago.ascasibar@uam.es>
Date: Fri, 4 Oct 2024 17:52:21 +0200
Subject: [PATCH] attribute name  and normalisation in

---
 src/pst/models.py                 |  6 +++---
 tutorials/models/analytical.ipynb | 30 +++++++++++++++++-------------
 2 files changed, 20 insertions(+), 16 deletions(-)

diff --git a/src/pst/models.py b/src/pst/models.py
index 6ddecea..b2919c2 100644
--- a/src/pst/models.py
+++ b/src/pst/models.py
@@ -37,7 +37,7 @@ def alpha_powerlaw(self, value):
 
     def ism_metallicity(self, times):
         m = self.stellar_mass_formed(times)
-        return self.ism_metallicity_today * np.power(m / m[-1], self.alpha_powerlaw)
+        return self.ism_metallicity_today * np.power(m / self.mass_today, self.alpha_powerlaw)
 
 def sfh_quenching_decorator(stellar_mass_formed):
     """A decorator for including a quenching event in a given SFH."""
@@ -403,7 +403,7 @@ def __init__(self, **kwargs):
         self.mass_norm = 1
         mtoday = self.stellar_mass_formed(self.today)
         self.mass_norm = self.mass_today / mtoday
-        self.metallicity_today =  kwargs['metallicity_today']
+        self.ism_metallicity_today =  kwargs['ism_metallicity_today']
 
     @u.quantity_input
     def stellar_mass_formed(self, times: u.Quantity):
@@ -414,7 +414,7 @@ def stellar_mass_formed(self, times: u.Quantity):
 
     @u.quantity_input
     def ism_metallicity(self, time : u.Gyr):
-        return np.full(time.size, fill_value=self.metallicity_today)
+        return np.full(time.size, fill_value=self.ism_metallicity_today)
 
 
 class LogNormalZPowerLawCEM(MassPropMetallicityMixin, LogNormalCEM):
diff --git a/tutorials/models/analytical.ipynb b/tutorials/models/analytical.ipynb
index 57247ac..b9fe386 100644
--- a/tutorials/models/analytical.ipynb
+++ b/tutorials/models/analytical.ipynb
@@ -115,7 +115,7 @@
    "source": [
     "today = 13 << u.Gyr\n",
     "mass_today = 1e11 << u.Msun\n",
-    "metallicity_today = 0.02 << u.dimensionless_unscaled\n",
+    "ism_metallicity_today = 0.02 << u.dimensionless_unscaled\n",
     "alpha_powerlaw = 1.5\n",
     "t0 = 20.0 << u.Gyr\n",
     "scale = 2.0\n",
@@ -138,11 +138,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#model = models.LogNormalQuenchedCEM(\n",
-    "model = models.LogNormalCEM(\n",
+    "model = models.LogNormalQuenchedCEM(\n",
+    "#model = models.LogNormalCEM(\n",
     "    today=today,\n",
     "    mass_today=mass_today, t0=t0, scale=scale, quenching_time=quenching_time,\n",
-    "    metallicity_today=metallicity_today, alpha_powerlaw=alpha_powerlaw,\n",
+    "    ism_metallicity_today=ism_metallicity_today, alpha_powerlaw=alpha_powerlaw,\n",
     ")"
    ]
   },
@@ -228,6 +228,16 @@
     "ssp_weights = model.interpolate_ssp_masses(ssp, t_obs=today)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d91db156-54a7-4ef2-8a0f-9316bc8e2322",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model.ism_metallicity(2*u.Gyr)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -245,7 +255,9 @@
     "\n",
     "mappable = ax.pcolormesh(ssp.ages.to_value(\"Gyr\"), ssp.metallicities, ssp_weights.to_value(u.Msun), norm=LogNorm())\n",
     "\n",
-    "plt.colorbar(mappable, ax=ax, label=r\"SSP mass (M$_\\odot$)\")"
+    "plt.colorbar(mappable, ax=ax, label=r\"SSP mass (M$_\\odot$)\")\n",
+    "\n",
+    "ax.plot(today-cosmic_time, model.ism_metallicity(cosmic_time), 'k-', alpha=.5)\n"
    ]
   },
   {
@@ -283,14 +295,6 @@
     "\n",
     "ax.plot(ssp.wavelength, sed)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "21a1fcf8-9faf-4339-b272-ccf636d3de49",
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {