can someone pls verify if this code is correct or note
import numpy as np import scipy.special as sp import matplotlib.pyplot as plt import pandas as pd from scipy.linalg import lstsq
def zernike_radial(n, m, r): """Compute the radial Zernike polynomial R_n^m(r).""" R = np.zeros_like(r) for k in range((n - abs(m)) // 2 + 1): c = (-1) ** k * sp.comb(n - k, k) * sp.comb(n - 2 * k, (n - abs(m)) // 2 - k) R += c * r ** (n - 2 * k) return R
def zernike(n, m, r, theta): """Compute the full Zernike polynomial.""" if m == 0: return zernike_radial(n, m, r) elif m > 0: return zernike_radial(n, m, r) * np.cos(m * theta) else: return zernike_radial(n, m, r) * np.sin(-m * theta)
def zernike_fit(x, y, z, max_order=4): """Fit Zernike polynomials to the surface deformation data.""" # Convert to polar coordinates r = np.sqrt(x2 + y2) theta = np.arctan2(y, x)
# Normalize r to be within [0, 1]
r /= np.max(r)
# Generate Zernike basis functions
terms = []
indices = []
for n in range(max_order + 1):
for m in range(-n, n + 1, 2): # Zernike indices
terms.append(zernike(n, m, r, theta))
indices.append((n, m))
Z = np.array(terms).T # Zernike basis matrix
# Solve for coefficients using least squares
coeffs, _, _, _ = lstsq(Z, z)
return coeffs, indices, Z
def load_excel(): """Prompt the user for an Excel file path and load X, Y, Z data.""" file_path = input("Enter the Excel file path (e.g., D:\aeg\ir telescope.xls): ")
try:
df = pd.read_excel(file_path, engine='xlrd' if file_path.endswith('.xls') else 'openpyxl')
x, y, z = df.iloc[:, 0], df.iloc[:, 1], df.iloc[:, 2] # First three columns
return np.array(x), np.array(y), np.array(z)
except Exception as e:
print(f"Error reading file: {e}")
return None, None, None
x, y, z = load_excel()
if x is not None and y is not None and z is not None: # Ask for Zernike polynomial order max_order = int(input("Enter the maximum Zernike order to fit (default is 4): ") or 4)
# Fit Zernike polynomials
coeffs, indices, Z = zernike_fit(x, y, z, max_order=max_order)
# Reconstruct the surface
z_fit = Z @ coeffs
# Calculate RMS values
rms_original = np.sqrt(np.mean(z**2))
rms_error = np.sqrt(np.mean((z - z_fit) ** 2))
rms_zernike_contribution = np.sqrt(np.sum(coeffs**2))
# Plot original and reconstructed surface
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
sc1 = ax[0].scatter(x, y, c=z, cmap='coolwarm')
ax[0].set_title("Original Surface Deformation")
fig.colorbar(sc1, ax=ax[0])
sc2 = ax[1].scatter(x, y, c=z_fit, cmap='coolwarm')
ax[1].set_title("Zernike Reconstructed Surface")
fig.colorbar(sc2, ax=ax[1])
plt.show()
# Print fitted coefficients
print("\nFitted Zernike Coefficients:")
for i, (n, m) in enumerate(indices):
print(f"Z({n},{m}) coefficient: {coeffs[i]:.6f}")
# Print RMS values
print("\nRMS Analysis:")
print(f"Original Surface RMS: {rms_original:.6e}")
print(f"Reconstruction RMS Error: {rms_error:.6e}")
print(f"Total RMS Contribution from Zernike Terms: {rms_zernike_contribution:.6e}")