Merge pull request #96 from mdecleir/plot_fit

Plot fitted models in A(lambda)/A(V)
karllark · Dec 16, 2021 · 393405f · 393405f
2 parents eff2378 + 1678fd4
commit 393405f
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 6 deletions.
diff --git a/measure_extinction/merge_obsspec.py b/measure_extinction/merge_obsspec.py
@@ -155,7 +155,9 @@ def merge_stis_obsspec(obstables, waveregion="UV", output_resolution=1000):
         cfluxes = ctable["FLUX"]
         cnpts = ctable["NPTS"].data
         for k in range(n_waves):
-            gvals = (cwaves >= full_wave_min[k]) & (cwaves < full_wave_max[k]) & (cnpts > 0)
+            gvals = (
+                (cwaves >= full_wave_min[k]) & (cwaves < full_wave_max[k]) & (cnpts > 0)
+            )
             if np.sum(gvals) > 0:
                 weights = 1.0 / np.square(cuncs[gvals].value)
                 full_flux[k] += np.sum(weights * cfluxes[gvals].value)
@@ -263,6 +265,7 @@ def merge_spex_obsspec(obstable, mask=[], output_resolution=2000):
     fluxes = obstable["FLUX"].data
     uncs = obstable["ERROR"].data
     npts = np.full((len(obstable["FLUX"])), 1.0)
+
     # take out data points that were "flagged" as bad by SpeXtool (i.e. FLAG is not zero)
     npts[obstable["FLAG"] != 0] = 0
     # take out data points with NaN fluxes

diff --git a/measure_extinction/plotting/plot_ext.py b/measure_extinction/plotting/plot_ext.py
@@ -215,7 +215,7 @@ def plot_extmodels(extdata, alax=False):
         plt.legend(bbox_to_anchor=(0.99, 0.9))
 
 
-def plot_fitmodel(extdata, yoffset=0, res=False):
+def plot_fitmodel(extdata, alax=False, yoffset=0, res=False):
     """
     Overplot a fitted model if available
 
@@ -224,6 +224,9 @@ def plot_fitmodel(extdata, yoffset=0, res=False):
     extdata : ExtData
         Extinction data under consideration
 
+    alax : boolean [default=False]
+        Whether or not to plot A(lambda)/A(X) instead of E(lambda-X)
+
     yoffset : float [default=0]
         Offset of the corresponding extinction curve (in order to match the model to the curve)
 
@@ -248,12 +251,19 @@ def plot_fitmodel(extdata, yoffset=0, res=False):
                 extdata.model["params"][0].value,
                 extdata.model["params"][2].value,
             )
-
         else:
             labeltxt = "fitted model"
+
+        # obtain the model extinctions
+        mod_ext = extdata.model["exts"]
+
+        # if the plot needs to be in A(lambda)/A(V), the model extinctions need to be converted to match the data
+        if alax:
+            mod_ext = (mod_ext / extdata.columns["AV"][0]) + 1
+
         plt.plot(
             extdata.model["waves"],
-            extdata.model["exts"] + yoffset,
+            mod_ext + yoffset,
             "-",
             lw=3,
             color="crimson",
@@ -522,7 +532,7 @@ def plot_multi_extinction(
 
         # overplot a fitted model if requested
         if fitmodel:
-            plot_fitmodel(extdata, yoffset=yoffset)
+            plot_fitmodel(extdata, alax=alax, yoffset=yoffset)
 
     # overplot Milky Way extinction curve models if requested
     if extmodels:
@@ -673,7 +683,7 @@ def plot_extinction(
 
     # overplot a fitted model if requested
     if fitmodel:
-        plot_fitmodel(extdata, res=True)
+        plot_fitmodel(extdata, alax=alax, res=True)
 
     # plot HI-lines if requested
     if HI_lines: