photometricoffsetfigures.ipynb [Emilys Galactic Gits]

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Analysis of Photometric Offsets in Astronomical Data using Python ======
=== PhotometricOffsetFigures.ipynb ===
This Jupyter notebook contains a detailed analysis focused on comparing photometric data from the HUGS and ACS surveys, specifically targeting the NGC 2808 cluster. It involves data cleaning, sampling, and visualization using libraries like pandas, matplotlib, and numpy. Furthermore, it applies fiducial line measurement techniques to discern photometric offsets and undertakes isochrone fitting with Bolometric Corrections to study stellar populations. The notebook concludes with a unique assessment of helium mass fraction effects on turn-off points in color-magnitude diagrams, potentially offering insights into stellar evolution.


```python
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from fidanka.fiducial.fiducial import measure_fiducial_lines
from fidanka.fiducial.utils import bin_color_mag_density
import pickle as pkl
from matplotlib.lines import Line2D
from fidanka.fiducial.fiducial import hull_density
from scipy.interpolate import splrep, BSpline
from scipy.interpolate import interp1d
from fidanka.isochrone.MIST import read_iso
from fidanka.isochrone.isochrone import interp_isochrone_age
from fidanka.bolometric.bctab import BolometricCorrector
import re
from scipy.optimize import curve_fit
from fidanka.isochrone.MIST import read_iso
from fidanka.isochrone.isochrone import interp_isochrone_age
```


```python
ACSRoot = "../photometry/hlsp_acsggct_hst_acs-wfc_ngc2808_r.rdviq.cal.adj.zpt"
HUGSRoot = "../photometry/HUGS/ngc2808/photometry.pkl"
```


```python
with open(HUGSRoot, 'rb') as f:
    HUGSdata = pkl.load(f)
```


```python
ACSdata = pd.read_csv(ACSRoot, sep=r'\s+')
ACSdata = ACSdata[ACSdata['Vvega'] - ACSdata['Ivega'] > 0.5]
```


```python
n = 20000
ACSSub = ACSdata.sample(n)
ACSmag = ACSSub["Vvega"]
ACScolor = ACSSub["Vvega"] - ACSSub["Ivega"]

HUGSSub = HUGSdata[1].sample(n)
HUGSmag = HUGSSub["F606W"]
HUGScolor = HUGSSub["F606W"] - HUGSSub["F814W"]
```


```python
ACSdensity = hull_density(ACScolor, ACSmag)
HUGSdensity = hull_density(HUGScolor, HUGSmag)
```


```python
with plt.style.context(pubStyle):
    fig, ax = plt.subplots(1,1,figsize=(10,7))
    ax.scatter(HUGScolor, HUGSmag,s=1,c=HUGSdensity, cmap="Reds")
    ax.scatter(ACScolor, ACSmag, s=1, c=ACSdensity, cmap='Blues', alpha=0.25)
    ax.set_xlim(0,2)
    ax.set_ylim(16, 25)
    ax.invert_yaxis()
    
    ax.set_xlabel("F606W - F814W", fontsize=27)
    ax.set_ylabel("F606W", fontsize=27)
    
    
    legend_elements = [
        Line2D([0], [0], marker='o', color='w', label='HUGS', markerfacecolor='r', markersize=15),
        Line2D([0], [0], marker='o', color='w', label='ACS', markerfacecolor='b', markersize=15),
    ]
    ax.legend(handles = legend_elements, frameon=False, fontsize=27)
    fig.savefig("Figures/photometricOffset.pdf")
```


    
![png](media/PhotometricOffsetFigures.ipynb/output_6_0.png)
    



```python
ACSFiducial = measure_fiducial_lines(
    ACSSub['Vvega'].values,
    ACSSub['Ivega'].values,
    ACSSub['err'].values,
    ACSSub['err.2'].values,
    binSize='uniformCS',
    targetStat=100,
    mcruns=5,
    convexHullPoints=50,
    nPops=1)
```

    100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:58<00:00, 11.61s/it]



```python
HUGSFiducial = measure_fiducial_lines(
    HUGSSub['F606W'].values,
    HUGSSub['F814W'].values,
    HUGSSub['F606W_RMS'].values,
    HUGSSub['F814W_RMS'].values,
    binSize='uniformCS',
    targetStat=100,
    mcruns=5,
    convexHullPoints=50,
    nPops=1)
```

    100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:52<00:00, 10.48s/it]



```python
from fidanka.misc.utils import measusre_perpendicular_distance
```


```python
with plt.style.context(pubStyle):
    fig, axs = plt.subplots(1,2,figsize=(20,7))
    ax = axs[0]
    ax.scatter(HUGScolor, HUGSmag,s=1,c=HUGSdensity, cmap="Reds")
    ax.scatter(ACScolor, ACSmag, s=1,c=ACSdensity, cmap='Blues', alpha=0.1)
    ax.set_xlim(0,2)
    ax.set_ylim(17, 22)
    ax.invert_yaxis()
    
    ACSSmoothTCK = splrep(ACSFiducial[0].mean[1], ACSFiducial[0].mean[0], s=0.001)
    ax.plot(BSpline(*ACSSmoothTCK)(ACSFiducial[0].mean[1]), ACSFiducial[0].mean[1] , color='g')
    HUGSmoothTCK = splrep(HUGSFiducial[0].mean[1], HUGSFiducial[0].mean[0], s=0.001)
    ax.plot(BSpline(*HUGSmoothTCK)(HUGSFiducial[0].mean[1]), HUGSFiducial[0].mean[1] , color='orange')
    ax.set_xlim(0.6, 1)
    
    ax.set_xlabel("F606W - F814W", fontsize=27)
    ax.set_ylabel("F606W", fontsize=27)
    
    
    legend_elements = [
        Line2D([0], [0], marker='o', color='w', label='HUGS', markerfacecolor='r', markersize=10),
        Line2D([0], [0], marker='o', color='w', label='ACS', markerfacecolor='b', markersize=10),
    ]
    ax.legend(handles = legend_elements, frameon=False, fontsize=27)
    
    
    ax = axs[1]
    ax.set_xlabel("F606W", fontsize=27)
    ax.set_ylabel("HUGS$_{color}$ - ACS$_{color}$", fontsize=27)
    
    domain = np.linspace(19.5, 22, 1000)
    dist = measusre_perpendicular_distance(BSpline(*HUGSmoothTCK), BSpline(*ACSSmoothTCK), domain)
    ax.plot(domain, dist, color='black')
    ax.axhline(np.mean(dist), color='black', linestyle='dashed')
    # ax.axhline(np.median(dist))
    fig.savefig("Figures/photometricOffset.pdf")
```


    
![png](media/PhotometricOffsetFigures.ipynb/output_10_0.png)
    



```python
plt.plot(np.arange(17, 22, 0.1),measusre_perpendicular_distance(BSpline(*HUGSmoothTCK), BSpline(*ACSSmoothTCK), np.arange(17, 22, 0.1)), 'o') 
```




    [<matplotlib.lines.Line2D at 0x7fa884d2c5b0>]




    
![png](media/PhotometricOffsetFigures.ipynb/output_11_1.png)
    



```python
bc = BolometricCorrector("WFC3", 1.13)
```


```python
rootPath = "../outputs.denseAlpha.fixedLowmass/"
```


```python
import pathlib
```


```python
paths = list()
for file in pathlib.Path(rootPath).rglob("alpha-1.901/isochrones.txt"):
    paths.append(file)
```


```python
isos = list()
testIso = read_iso(paths[0])
header = testIso[list(testIso.keys())[0]].columns
for path in paths:
    isos.append(pd.DataFrame(interp_isochrone_age(read_iso(path), 12.3), columns=header))
```


```python
isos[0]
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>EEP</th>
      <th>log10_isochrone_age_yr</th>
      <th>initial_mass</th>
      <th>star_mass</th>
      <th>star_mdot</th>
      <th>he_core_mass</th>
      <th>c_core_mass</th>
      <th>log_L</th>
      <th>log_LH</th>
      <th>log_LHe</th>
      <th>...</th>
      <th>surface_c12</th>
      <th>surface_o16</th>
      <th>log_center_T</th>
      <th>log_center_Rho</th>
      <th>center_gamma</th>
      <th>center_h1</th>
      <th>center_he4</th>
      <th>center_c12</th>
      <th>phase</th>
      <th>age</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>286.0</td>
      <td>10.089905</td>
      <td>0.235747</td>
      <td>0.235747</td>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.0</td>
      <td>-1.989956</td>
      <td>-1.989276</td>
      <td>-10.0</td>
      <td>...</td>
      <td>0.000061</td>
      <td>0.000820</td>
      <td>6.849942</td>
      <td>2.109066</td>
      <td>0.0</td>
      <td>0.688852</td>
      <td>0.309676</td>
      <td>0.000061</td>
      <td>0.0</td>
      <td>1.230187e+10</td>
    </tr>
    <tr>
      <th>1</th>
      <td>287.0</td>
      <td>10.089905</td>
      <td>0.260267</td>
      <td>0.260267</td>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.0</td>
      <td>-1.904474</td>
      <td>-1.903862</td>
      <td>-10.0</td>
      <td>...</td>
      <td>0.000061</td>
      <td>0.000820</td>
      <td>6.863070</td>
      <td>2.051913</td>
      <td>0.0</td>
      <td>0.688946</td>
      <td>0.309582</td>
      <td>0.000061</td>
      <td>0.0</td>
      <td>1.230187e+10</td>
    </tr>
    <tr>
      <th>2</th>
      <td>288.0</td>
      <td>10.089905</td>
      <td>0.262639</td>
      <td>0.262639</td>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.0</td>
      <td>-1.897316</td>
      <td>-1.896708</td>
      <td>-10.0</td>
      <td>...</td>
      <td>0.000061</td>
      <td>0.000820</td>
      <td>6.864233</td>
      <td>2.047220</td>
      <td>0.0</td>
      <td>0.688909</td>
      <td>0.309619</td>
      <td>0.000061</td>
      <td>0.0</td>
      <td>1.230187e+10</td>
    </tr>
    <tr>
      <th>3</th>
      <td>289.0</td>
      <td>10.089905</td>
      <td>0.265015</td>
      <td>0.265015</td>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.0</td>
      <td>-1.890155</td>
      <td>-1.889553</td>
      <td>-10.0</td>
      <td>...</td>
      <td>0.000061</td>
      <td>0.000820</td>
      <td>6.865400</td>
      <td>2.042532</td>
      <td>0.0</td>
      <td>0.688879</td>
      <td>0.309649</td>
      <td>0.000061</td>
      <td>0.0</td>
      <td>1.230187e+10</td>
    </tr>
    <tr>
      <th>4</th>
      <td>290.0</td>
      <td>10.089905</td>
      <td>0.267398</td>
      <td>0.267398</td>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.0</td>
      <td>-1.882988</td>
      <td>-1.882392</td>
      <td>-10.0</td>
      <td>...</td>
      <td>0.000061</td>
      <td>0.000820</td>
      <td>6.866570</td>
      <td>2.037846</td>
      <td>0.0</td>
      <td>0.688857</td>
      <td>0.309671</td>
      <td>0.000061</td>
      <td>0.0</td>
      <td>1.230187e+10</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>187</th>
      <td>601.0</td>
      <td>10.089905</td>
      <td>0.763934</td>
      <td>0.412178</td>
      <td>0.0</td>
      <td>0.332015</td>
      <td>0.0</td>
      <td>2.349763</td>
      <td>2.347601</td>
      <td>-10.0</td>
      <td>...</td>
      <td>0.000056</td>
      <td>0.000788</td>
      <td>7.712642</td>
      <td>5.599579</td>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.998341</td>
      <td>0.000004</td>
      <td>2.0</td>
      <td>1.230187e+10</td>
    </tr>
    <tr>
      <th>188</th>
      <td>602.0</td>
      <td>10.089905</td>
      <td>0.763948</td>
      <td>0.410825</td>
      <td>0.0</td>
      <td>0.333765</td>
      <td>0.0</td>
      <td>2.364422</td>
      <td>2.362240</td>
      <td>-10.0</td>
      <td>...</td>
      <td>0.000056</td>
      <td>0.000788</td>
      <td>7.715019</td>
      <td>5.605118</td>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.998341</td>
      <td>0.000004</td>
      <td>2.0</td>
      <td>1.230187e+10</td>
    </tr>
    <tr>
      <th>189</th>
      <td>603.0</td>
      <td>10.089905</td>
      <td>0.763962</td>
      <td>0.409429</td>
      <td>0.0</td>
      <td>0.335503</td>
      <td>0.0</td>
      <td>2.379073</td>
      <td>2.376870</td>
      <td>-10.0</td>
      <td>...</td>
      <td>0.000056</td>
      <td>0.000788</td>
      <td>7.717450</td>
      <td>5.610761</td>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.998342</td>
      <td>0.000004</td>
      <td>2.0</td>
      <td>1.230187e+10</td>
    </tr>
    <tr>
      <th>190</th>
      <td>604.0</td>
      <td>10.089905</td>
      <td>0.763978</td>
      <td>0.407971</td>
      <td>0.0</td>
      <td>0.337333</td>
      <td>0.0</td>
      <td>2.393712</td>
      <td>2.391482</td>
      <td>-10.0</td>
      <td>...</td>
      <td>0.000056</td>
      <td>0.000788</td>
      <td>7.719974</td>
      <td>5.616615</td>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.998342</td>
      <td>0.000004</td>
      <td>2.0</td>
      <td>1.230187e+10</td>
    </tr>
    <tr>
      <th>191</th>
      <td>605.0</td>
      <td>10.089905</td>
      <td>0.764000</td>
      <td>0.406203</td>
      <td>0.0</td>
      <td>0.339454</td>
      <td>0.0</td>
      <td>2.408234</td>
      <td>2.405955</td>
      <td>-10.0</td>
      <td>...</td>
      <td>0.000056</td>
      <td>0.000788</td>
      <td>7.722938</td>
      <td>5.623576</td>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.998342</td>
      <td>0.000004</td>
      <td>3.0</td>
      <td>1.230187e+10</td>
    </tr>
  </tbody>
</table>
<p>192 rows × 26 columns</p>
</div>




```python
mags = list()
for iso in isos:
    mags.append(bc.apparent_mags(
        10**iso['log_Teff'],
        iso['log_g'],
        iso['log_L'],
        filters=("WFC3_UVIS_F606W", "WFC3_UVIS_F814W")
    )
               )
```


```python
turnOffs = list()
for mag, path in zip(mags, paths):
    color = mag["WFC3_UVIS_F606W"]-mag["WFC3_UVIS_F814W"]
    idStr = re.findall('Pop[A|E]\+0\.\d+', str(path))[0]
    pop = idStr[3]
    Y = float(idStr[5:])
    turnOffColor = min(mag["WFC3_UVIS_F606W"] - mag["WFC3_UVIS_F814W"])
    turnOffMag = mag[color==turnOffColor]["WFC3_UVIS_F606W"]
    if turnOffColor < 0.48:
        turnOffs.append((turnOffColor, turnOffMag.values[0], pop, Y))
    
meanTurnOffColorA = sum([x[0] for x in turnOffs if x[2] == 'A'])/(len(turnOffs)/2)
meanTurnOffMagA = sum([x[1] for x in turnOffs if x[2] == 'A'])/(len(turnOffs)/2)
meanTurnOffColorE = sum([x[0] for x in turnOffs if x[2] == 'E'])/(len(turnOffs)/2)
meanTurnOffMagE = sum([x[1] for x in turnOffs if x[2] == 'E'])/(len(turnOffs)/2)
```


```python
hevolve = list()
for mag, path, turnOff in zip(mags, paths, turnOffs):
    formater = {'A': 'ko', 'E': 'kx'}
    if turnOff[2] == 'A':
        dist = np.sqrt((turnOff[0]-turnOffs[0][0])**2 + (turnOff[1]-turnOffs[0][1])**2)
        hevolve.append((turnOff[3]-turnOffs[0][3], dist))
hevolve = np.array(hevolve)
```


```python
with plt.style.context(pubStyle):
    fig, ax = plt.subplots(1,1,figsize=(10,7))
    for h in hevolve:
        ax.plot(h[0], h[1], 'ko')
        
    line = lambda x, m, b:  m*x + b
    fit, covar = curve_fit(line, hevolve[:, 0], hevolve[:, 1])
    X = np.linspace(hevolve[:, 0].min(), hevolve[:, 0].max())
    ax.plot(X, line(X, *fit), 'k--')
    ax.set_xlabel("$\Delta$ Helium Mass Fraction", fontsize=20)
    ax.set_ylabel("Turn off offset from mean [mag]", fontsize=20)
    ax.annotate(f'$O(Y)={fit[0]:0.2f}Y+{fit[1]:0.2f}$', xy=(-0.04, 0.30), color='black', fontsize=20)

    fig.savefig("Figures/HeliumMeanOffset.pdf")
```


    
![png](media/PhotometricOffsetFigures.ipynb/output_21_0.png)
    



```python
with open("BestFit.pkl", 'rb') as f:
    mags = pkl.load(f)
```


```python
def correct_iso_set(path, age, mu, Av):
    print(os.path.exists(path))
    
    isoAtAge = interp_isochrone_age(iso, age)

    isoAMags = bc.apparent_mags(
        10**isoAtAge[:, 10],
        isoAtAge[:, 12],
        isoAtAge[:, 7],
        mu = mu,
        Av = Av,
        filters = ("WFC3_UVIS_F275W", "WFC3_UVIS_F606W", "WFC3_UVIS_F814W")
    )
    return isoAMags
```


```python
testPath = "/mnt/Astronomy/GraduateSchool/Thesis/GCConsistency/NGC2808/justIso.denseAlpha.fixedLowmass/PopA+0.39/alpha-1.750/isochrones.txt"
```


```python
isoMags = correct_iso_set(testPath, 17, 1.491e+01, 5.414e-01)
```

    True



```python
fig, ax = plt.subplots(1,1,figsize=(10,10))
ax.scatter(HUGScolor, HUGSmag,s=1,c=HUGSdensity, cmap="Reds")
ax.scatter(ACScolor, ACSmag, s=1,c=ACSdensity, cmap='Blues', alpha=0.1)
ax.set_ylim(17, 22)
ax.invert_yaxis()
ACSSmoothTCK = splrep(ACSFiducial[0].mean[1], ACSFiducial[0].mean[0], s=0.001)
ax.plot(BSpline(*ACSSmoothTCK)(ACSFiducial[0].mean[1]), ACSFiducial[0].mean[1] , color='C0')
HUGSmoothTCK = splrep(HUGSFiducial[0].mean[1], HUGSFiducial[0].mean[0], s=0.001)
ax.plot(BSpline(*HUGSmoothTCK)(HUGSFiducial[0].mean[1]), HUGSFiducial[0].mean[1] , color='C1')
ax.set_xlim(0.6, 1)

ax.set_xlabel("F606W - F814W", fontsize=27)
ax.set_ylabel("F606W", fontsize=27)

Color = isoMags['WFC3_UVIS_F606W'].values - isoMags["WFC3_UVIS_F814W"].values
Mag = isoMags['WFC3_UVIS_F606W'].values
turnOffColor = min(Color)
turnOffMag = Mag[Color==turnOffColor][0]
print(turnOffMag, turnOffColor)

ax.plot(Color, Mag, color='C2', linewidth=10)
```

    19.73960461108486 0.652611121597868





    [<matplotlib.lines.Line2D at 0x7fa87ae821c0>]




    
![png](media/PhotometricOffsetFigures.ipynb/output_26_2.png)
    



```python

```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>WFC3_UVIS_F275W</th>
      <th>WFC3_UVIS_F606W</th>
      <th>WFC3_UVIS_F814W</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>33.753726</td>
      <td>26.699337</td>
      <td>24.461482</td>
    </tr>
    <tr>
      <th>1</th>
      <td>33.737848</td>
      <td>26.683081</td>
      <td>24.448594</td>
    </tr>
    <tr>
      <th>2</th>
      <td>33.722102</td>
      <td>26.666965</td>
      <td>24.435789</td>
    </tr>
    <tr>
      <th>3</th>
      <td>33.706476</td>
      <td>26.650976</td>
      <td>24.423060</td>
    </tr>
    <tr>
      <th>4</th>
      <td>33.690944</td>
      <td>26.635087</td>
      <td>24.410393</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>313</th>
      <td>23.471173</td>
      <td>14.562014</td>
      <td>13.486965</td>
    </tr>
    <tr>
      <th>314</th>
      <td>23.445802</td>
      <td>14.533347</td>
      <td>13.452625</td>
    </tr>
    <tr>
      <th>315</th>
      <td>23.420594</td>
      <td>14.504854</td>
      <td>13.418372</td>
    </tr>
    <tr>
      <th>316</th>
      <td>23.395397</td>
      <td>14.476370</td>
      <td>13.384128</td>
    </tr>
    <tr>
      <th>317</th>
      <td>23.370260</td>
      <td>14.447949</td>
      <td>13.349919</td>
    </tr>
  </tbody>
</table>
<p>318 rows × 3 columns</p>
</div>




```python

```