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
Define Zernike polynomial functions
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
Function to load data from an Excel file
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
Get user input for the Excel file path
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}")
