I'm unable to figure out what to do, I can't seem to identify the minima's error here. I have another dataset where I get my variables from.
# Adjusting the code to correctly identify the minimum of the polynomial fit.
import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial
# Load the data from the file
data = np.genfromtxt('C:/Users/me/Desktop/alex 2/V0520_cas/23-10-11/results_hjd.diff', dtype=None, encoding=None)
hjd, mag, mag_err = data[:, 0], data[:, 1], data[:, 2]
# Replace zero or negative errors with the median of the positive errors
mag_err[mag_err <= 0] = np.median(mag_err[mag_err > 0])
# Define the window size and get the subset of data around the 342nd data point
index_342 = 342
window_size = 30 # The number of data points to include from each side
# Select the data around the index_342
subset_slice = slice(max(index_342 - window_size, 0), min(index_342 + window_size + 1, len(hjd)))
hjd_subset, mag_subset, mag_err_subset = hjd[subset_slice], mag[subset_slice], mag_err[subset_slice]
# Fit a polynomial of degree 4 to the data
coeffs = np.polyfit(hjd_subset, mag_subset, 4, w=1/mag_err_subset)
# Create a polynomial object from the coefficients
p = np.poly1d(coeffs)
# Generate points for the fitted curve
hjd_fit = np.linspace(hjd_subset.min(), hjd_subset.max(), 500)
mag_fit = p(hjd_fit)
# Find the minimum point on the curve which corresponds to the maximum light
# We'll use np.polyder to find the derivative of the polynomial and np.roots to find its roots.
derivative = np.polyder(p)
# Roots of the derivative (critical points)
critical_points = np.roots(derivative)
# Filter roots to include only those within the range of our subset data
critical_points = critical_points[(critical_points > hjd_subset.min()) & (critical_points < hjd_subset.max())]
# Evaluate the polynomial at the critical points to find the minimum
critical_mags = p(critical_points)
min_index = np.argmin(critical_mags)
min_hjd = critical_points[min_index]
max_mag = critical_mags[min_index]
# Plotting the data and the fitted polynomial curve
plt.figure(figsize=(12, 6))
plt.errorbar(hjd_subset, mag_subset, yerr=mag_err_subset, fmt='o', label='Data Subset', markersize=5)
plt.plot(hjd_fit, mag_fit, label='Polynomial Fit', color='red')
plt.axvline(x=min_hjd, color='green', linestyle='--', label=f'Minimum Light at HJD {min_hjd:.5f}')
plt.xlabel('Heliocentric Julian Date (HJD)')
plt.ylabel('Magnitude')
plt.title('Light Curve with Polynomial Fit')
plt.gca().invert_yaxis() # Invert y-axis for astronomical convention
plt.legend()
plt.show()
# Output the time of minimum light (maximum light)
print(f"Time of Max Light (HJD): {min_hjd:.5f}")
print(f"Magnitude at Max Light: {max_mag:.5f}")
I tried to estimate error from stdev of residuals as follows:
# Estimate the error of the maximum light from the standard deviation of the residuals
residuals = mag_subset - p_fit(hjd_subset)
residual_std = np.std(residuals)
# Plotting the data and the fitted polynomial curve
plt.figure(figsize=(12, 6))
plt.errorbar(hjd_subset, mag_subset, yerr=mag_err_subset, fmt='o', label='Data Subset', markersize=5)
plt.plot(hjd_fit, mag_fit, label='Polynomial Fit', color='red')
plt.axvline(x=min_hjd, color='green', linestyle='--', label=f'Maximum Light at HJD {min_hjd:.5f}')
plt.xlabel('Heliocentric Julian Date (HJD)')
plt.ylabel('Magnitude')
plt.title('Light Curve with Polynomial Fit')
plt.gca().invert_yaxis() # Invert y-axis for astronomical convention
plt.legend()
plt.show()
# Output the time of maximum light and its error
print(f"Time of Maximum Light (HJD): {min_hjd:.5f}")
print(f"Error in Time of Maximum Light (Estimated): {residual_std:.5f}")
print(f"Magnitude at Maximum Light: {max_mag:.5f}")
enter image description here I get this from my error calculations clearly wrong
