How to centre observational data onto predicted isochrone data in Python

47 views Asked by At

Edit: The observational data can be found here. The isochrone or predicted data can be found here.

So I'm trying to best fit observational astronomy data to a predicted 'patttern' callen an isochrone. The isochrone isn't a function and merely a prediction from calculations. This is shown in the image below:

Obverved data vs. predicted isochrone

The goal is to make x and y adjustments to the predicted data to make a best fit. Not all parts of the predicted data is necessarily represented in the observed data, but the particular part in the data roughly between 2 < y < 6 has a nice curve which is present in all isochrones, and at the very least that's the part I'd like to fit.

The issue is that since the predicted isochrone isn't a function, I can't use scipy.curvefit.

I have attempted to fit a polynomial to the area of interest in the isochrone, and used that with scipy.curvefit to get some horizontal displacement.

poly_coeffs = np.polyfit(iso_g, iso_bprp, deg=40)
poly = np.poly1d(poly_coeffs)

### Plot these sideways
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot()

ax.plot(iso_bprp, iso_g)
ax.plot(poly(iso_g), iso_g)

ax.set_xlabel(r'Colour Index [$Bp - Rp$]', fontsize=20)
ax.set_ylabel('G Mag', fontsize=20)

ax.invert_yaxis()
ax.tick_params(axis='both', which='major', labelsize=15)
plt.show()

Fitting the isochrone to a polynomial

I was able to make a shift in this example:

Horizontal shift

However, I can't create a vertical shift when try to introduce that parameter:

def fit(x, xshift, yshift):
    return poly(x + xshift) + yshift

parameters, covariance = curve_fit(fit, data_g, data_bprp)

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot()

ax.plot(poly(iso_g), iso_g, label='Original Curve')
ax.scatter(data_bprp, data_g, s=0.5)

ax.plot(fit(iso_g, parameters[0], parameters[1]), iso_g, label='XY Shifted Curve')

ax.set_xlabel(r'Colour Index [$Bp - Rp$]', fontsize=20)
ax.set_ylabel('G Mag', fontsize=20)

ax.invert_yaxis()
ax.tick_params(axis='both', which='major', labelsize=15)
plt.legend()
plt.show()

Attempted XY shift

I'm sure there should be a better way to acheive this goal, but I am unable to find it. Any insight is appreciated. Thanks!

0

There are 0 answers