import numpy as np
import plotly.graph_objects as go
from sklearn.linear_model import LinearRegression

# -------------------------------
# Synthetic dataset (replace with experimental data)
# -------------------------------
N = 150
frequencies = np.linspace(100, 1000, N)
X = np.vstack([
    frequencies,
    frequencies**2,
    np.sin(frequencies/100),
    np.cos(frequencies/200),
    np.power(1.618, frequencies/500),
    np.sqrt(frequencies),
    np.log1p(frequencies)
]).T

Y = np.vstack([
    3*np.sin(frequencies/200)+0.1*np.random.randn(N),
    6*np.sin(frequencies/400)+0.1*np.random.randn(N),
    2*np.cos(frequencies/300)+0.05*np.random.randn(N),
    3*np.cos(frequencies/350)+0.05*np.random.randn(N)
]).T

# Fit linear model (T matrix)
model = LinearRegression().fit(X, Y)

# Function to predict cymatic params from frequency
def predict_params(f):
    features = np.array([
        f,
        f**2,
        np.sin(f/100),
        np.cos(f/200),
        np.power(1.618, f/500),
        np.sqrt(f),
        np.log1p(f)
    ]).reshape(1,-1)
    pred = model.predict(features)[0]
    return dict(alpha=abs(pred[0]), beta=abs(pred[1]), eta=abs(pred[2]), zeta=abs(pred[3]))

# -------------------------------
# Generate Smith-chart-style polar points
# -------------------------------
freq_sweep = np.linspace(100, 1000, 100)
radii = []
thetas = []

for f in freq_sweep:
    p = predict_params(f)
    # Radius = amplitude scaling (η)
    r = p['eta']
    # Angle = mode combination (α + β + ζ) scaled to 2pi
    theta = ((p['alpha'] + p['beta'] + p['zeta']) % 2) * np.pi
    radii.append(r)
    thetas.append(theta)

# Convert to Cartesian for Plotly scatterpolar
x = np.array(radii) * np.cos(thetas)
y = np.array(radii) * np.sin(thetas)

# Plot interactive Smith-like chart
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y, mode='lines+markers', line=dict(color='blue'), marker=dict(size=6)))
fig.update_layout(
    title='Cymatic Smith Chart (Polar Resonance Trajectory)',
    xaxis_title='Re (Amplitude cosθ)',
    yaxis_title='Im (Amplitude sinθ)',
    xaxis=dict(scaleanchor='y', scaleratio=1),
    yaxis=dict(constrain='domain'),
    showlegend=False
)
fig.show()
