Let’s start performing some data analysis with our beautiful JWST near-infrared catalogs. But first! let’s go over what color-magnitude diagrams (CMDs) are. CMDs are the life and blood for astronomers. If you are familiar with the Hertzsprung Russell (HR) diagram, CMDs are a subset of these plots.

They are actually fairly simple. You merely take magnitudes of stars from one catalog ($y$) and plot it relative to the difference in magnitudes between two different catalogs ($x$). An example for both F115W vs. (F115W - F212N) and F212N vs. (F115W - F212N) is shown below.

It is standardized to plot the (smaller $\lambda$ - larger $\lambda$) on the $x$-axis.

Now without getting into too much detail, different types
of stars occupy different regions of the CMD. The large dense cluster of stars
in the plots above showcase the *main-sequence* track of the CMD, and the
linear streak above the main-sequence cluster contains the *Red Clump (RC)*
stars – which are the stars I will be using in my research.

Let’s actually build the above CMDs for the F115W and F212N catalogs. As the
catalogs are of different size, with a vastly different # of stars, we first
have to determine which stars are *common* to both catalogs.

As the catalogs are close in wavelength and were taken at approximately the same time, we can assume these catalogs place the star’s $x$ and $y$ coordinates on the same coordinate system. Indeed, if we overplot both the F212N and F115W catalog centroid coordinates on top of each other, with an very zoomed in image in the background…

Here, the red squares indicate stars from the *F212N* catalog, and the green
squares indicate stars from the *F115W* catalog.

As there are `121349`

F212N
stars but only `53620`

F115W stars for this image, there are a lot more red
markers than green. However, if you look closely you can see for the very
luminous stars, the green and red markers exist on top of one another. Indeed,
the two catalogs are on roughly the same coordinate system. (*Note* for future work,
I do perform coordinate transformations on all my catalogs…this is
preliminary work*).

Finally, with all the introductions out of the way, we can start doing some research.

*Finding Common Stars Across F115W & F212N*

To determine the stars common to both F115w and F212n, I employ a nearest-neighbor KDTree algorithm.

No transformations are done and it is assumed that the two catalogs are already on the same coordinate system and magnitude system.

For two stars to be matched, they must be within a specified radius (`dr_tol`

) and
delta-magnitude (`dm_tol`

). For stars with more than 1 neighbor (within the tolerances),
if one is found that is the best match in both brightness and positional offsets
(closest in both), then the match is made. Otherwise,
their is a conflict and no match is returned for the star.

The entire code is shown on my GitHub, but I will lay out some here.

I first start by making coordinate pairs for each star in both catalogs:

```
if x1.shape != y1.shape:
raise ValueError('x1 and y1 do not match!')
if x2.shape != y2.shape:
raise ValueError('x2 and y2 do not match!')
# Setup coords1 pairs and coords 2 pairs
# this is equivalent to, but faster than just doing np.array([x1, y1])
coords1 = np.empty((x1.size, 2))
coords1[:, 0] = x1
coords1[:, 1] = y1
# this is equivalent to, but faster than just doing np.array([x1, y1])
coords2 = np.empty((x2.size, 2))
coords2[:, 0] = x2
coords2[:, 1] = y2
```

I then implement a KDTree using `scipy.spatial`

```
kdt = KDT(coords2, balanced_tree=False)
# This returns the number of neighbors within the specified
# radius. We will use this to find those stars that have no or one
# match and deal with them easily. The more complicated conflict
# cases will be dealt with afterward.
i2_match = kdt.query_ball_point(coords1, dr_tol)
Nmatch = np.array([len(idxs) for idxs in i2_match])
```

Afterwards, I loop through and handle all the different matches within `dr_tol`

and `dm_tol`

and place the indices of each star in NumPy Arrays `idxs1`

and `idxs2`

,
for the first and second catalogs respectively. I then deal with duplicates

```
# Deal with duplicates
duplicates = [item for item, count in list(Counter(idxs2).items()) if count > 1]
if verbose > 2:
print(( ' Found {0:d} duplicates out of {1:d} matches'.format(len(duplicates), len(dm))))
keep = np.ones(len(idxs1), dtype=bool)
for dd in range(len(duplicates)):
# Index into the idxs1, idxs2 array of this duplicate.
dups = np.where(idxs2 == duplicates[dd])[0]
# Assume the duplicates are confused first... see if we
# can resolve the confusion below.
keep[dups] = False
dm_dups = m1[idxs1[dups]] - m2[idxs2[dups]]
dr_dups = np.hypot(x1[idxs1[dups]] - x2[idxs2[dups]], y1[idxs1[dups]] - y2[idxs2[dups]])
dm_min = np.abs(dm_dups).argmin()
dr_min = np.abs(dr_dups).argmin()
# If there is a clearly preferred match (closest in distance and brightness), then
# keep it and dump the other duplicates.
if dm_min == dr_min:
keep[dups[dm_min]] = True
else:
if verbose:
print('confused, dropping')
# Clean up the duplicates
idxs1 = idxs1[keep]
idxs2 = idxs2[keep]
```

I use the above algorithm with the F115w and F212n catalogs, with a `dr_tol`

of
0.5 and a `dm_tol`

of 99 (I don’t care about `dm_tol`

right now).

And since the algorithm returns the indices of matched stars from their
respective catalogs, I can get the magnitudes for the matched stars by just
taking the values from the `m`

column when the index is in `idxs1`

or `idxs2`

.

```
x1 = df['x'] # x centroid position of catalog1 stars
y1 = df['y'] # y centroid position of catalog1 stars
m1 = df['m'] # vega magnitude of catalog1 stars
me1 = df['me'] # error in vega magnitude of catalog1 stars
x2 = df2['x']
y2 = df2['y']
m2 = df2['m']
me2 = df2['me']
m1_error = []
m2_error = []
m1_matched = []
m2_matched = []
m_difference = []
for i in idxs1:
m1_error += [me1[i]]
m1_matched += [m1[i]]
for i in idxs2:
m2_error += [me2[i]]
m2_matched += [m2[i]]
```

And since the CMDs are a plot of a magnitude vs mag-difference, I subtract the
`m1_matched`

and `m2_matched`

NumPy arrays to get the $x$-axis for all the
CMDs.

```
m_difference = np.subtract(m1_matched, m2_matched)
```

Now we have everything we need to generate the CMDs for our catalogs. I implement every step above for all the NRCB1 - NRCB4 catalogs for F115W and F212N.

```
N1_115_catalog = Table.read("catalogs/dr2/NRCB1_catalog115w.csv").to_pandas()
N1_212_catalog = Table.read("catalogs/dr2/NRCB1_catalog212n.csv").to_pandas()
N1_115_matched, N1_212_matched, N1_115_error, N1_212_error, N1_115_212_difference = matched(N1_115_catalog, N1_212_catalog, "NRCB1_F115W", "NRCB1_F212N", 0.5, 99)
N2_115_catalog = Table.read("catalogs/dr2/NRCB2_catalog115w.csv").to_pandas()
N2_212_catalog = Table.read("catalogs/dr2/NRCB2_catalog212n.csv").to_pandas()
N2_115_matched, N2_212_matched, N2_115_error, N2_212_error, N2_115_212_difference = matched(N2_115_catalog, N2_212_catalog, "NRCB2_F115W", "NRCB2_F212N", 0.5, 99)
N3_115_catalog = Table.read("catalogs/dr2/NRCB3_catalog115w.csv").to_pandas()
N3_212_catalog = Table.read("catalogs/dr2/NRCB3_catalog212n.csv").to_pandas()
N3_115_matched, N3_212_matched, N3_115_error, N3_212_error, N3_115_212_difference = matched(N3_115_catalog, N3_212_catalog, "NRCB3_F115W", "NRCB3_F212N", 0.5, 99)
N4_115_catalog = Table.read("catalogs/dr2/NRCB4_catalog115w.csv").to_pandas()
N4_212_catalog = Table.read("catalogs/dr2/NRCB4_catalog212n.csv").to_pandas()
N4_115_matched, N4_212_matched, N4_115_error, N4_212_error, N4_115_212_difference = matched(N4_115_catalog, N4_212_catalog, "NRCB4_F115W", "NRCB4_F212N", 0.5, 99)
```

Generating the CMDs

```
fig, axis = plt.subplots(4, 2, figsize = (20, 20))
fig.suptitle(' NRCB1-NRCB4 All Color-Magnitude Diagrams ', fontsize=15)
axis[0,0].invert_yaxis()
axis[0,1].invert_yaxis()
axis[1,0].invert_yaxis()
axis[1,1].invert_yaxis()
axis[2,0].invert_yaxis()
axis[2,1].invert_yaxis()
axis[3,0].invert_yaxis()
axis[3,1].invert_yaxis()
xy = np.vstack([N1_115_212_difference, N1_115_matched])
z = gaussian_kde(xy)(xy)
axis[0,0].scatter(N1_115_212_difference, N1_115_matched, c = z, s = 1)
axis[0,0].set_xlabel('NRCB1 F115W - F212N')
axis[0,0].set_ylabel('NRCB1 F115W')
axis[0,0].set_title('NRCB1 F115W vs. F115W - F212N')
xy = np.vstack([N1_115_212_difference, N1_212_matched])
z = gaussian_kde(xy)(xy)
axis[0,1].scatter(N1_115_212_difference, N1_212_matched, c = z, s = 1)
axis[0,1].set_xlabel('NRCB1 F115W - F212N')
axis[0,1].set_ylabel('NRCB1 F212N')
axis[0,1].set_title('NRCB1 F212N vs. F115W - F212N')
xy = np.vstack([N2_115_212_difference, N2_115_matched])
z = gaussian_kde(xy)(xy)
axis[1,0].scatter(N2_115_212_difference, N2_115_matched, c = z, s = 1)
axis[1,0].set_xlabel('NRCB2 F115W - F212N')
axis[1,0].set_ylabel('NRCB2 F115W')
axis[1,0].set_title('NRCB2 F115W vs. F115W - F212N')
xy = np.vstack([N2_115_212_difference, N2_212_matched])
z = gaussian_kde(xy)(xy)
axis[1,1].scatter(N2_115_212_difference, N2_212_matched, c = z, s = 1)
axis[1,1].set_xlabel('NRCB2 F115W - F212N')
axis[1,1].set_ylabel('NRCB2 F212N')
axis[1,1].set_title('NRCB2 F212N vs. F115W - F212N')
xy = np.vstack([N3_115_212_difference, N3_115_matched])
z = gaussian_kde(xy)(xy)
axis[2,0].scatter(N3_115_212_difference, N3_115_matched, c = z, s = 1)
axis[2,0].set_xlabel('NRCB3 F115W - F212N')
axis[2,0].set_ylabel('NRCB3 F115W')
axis[2,0].set_title('NRCB3 F115W vs. F115W - F212N')
xy = np.vstack([N3_115_212_difference, N3_212_matched])
z = gaussian_kde(xy)(xy)
axis[2,1].scatter(N3_115_212_difference, N3_212_matched, c = z, s = 1)
axis[2,1].set_title('NRCB3 F212N vs. F115W - F212N')
axis[2,1].set_xlabel('NRCB3 F115W - F212N')
axis[2,1].set_ylabel('NRCB3 F212N')
xy = np.vstack([N4_115_212_difference, N4_115_matched])
z = gaussian_kde(xy)(xy)
axis[3,0].scatter(N4_115_212_difference, N4_115_matched, c = z, s = 1)
axis[3,0].set_xlabel('NRCB4 F115W - F212N')
axis[3,0].set_ylabel('NRCB4 F115W')
axis[3,0].set_title('NRCB4 F115W vs. F115W - F212N')
xy = np.vstack([N4_115_212_difference, N4_212_matched])
z = gaussian_kde(xy)(xy)
axis[3,1].scatter(N4_115_212_difference, N4_212_matched, c = z, s = 1)
axis[3,1].set_xlabel('NRCB4 F115W - F212N')
axis[3,1].set_ylabel('NRCB4 F212N')
axis[3,1].set_title('NRCB4 F212N vs. F115W - F212N')
fig.tight_layout()
```

Beautiful! These CMDs have a clear main-sequence track and a very prominent RC bar with a definite slope. The NRCB1 (top row) CMDs match exactly to the CMDs we saw originally. The region NRCB3 is of particular interest since its the region containing Sagittarius A*.

We will now work on seeing how well past extinction laws match the slope of the RC Bar by utilizing theoretical synthetic isochrones for the Galactic Center. Eventually, we’ll derive our own, newer, more accurate extinction laws.