"""Independent polar-chart science reference; not a second browser starspot Lab."""
import numpy as np
from numpy.polynomial.legendre import leggauss


def disk_quadrature(radial_order=64, azimuth_order=240):
    nodes, weights = leggauss(radial_order)
    mu = (nodes + 1) / 2
    phi = (np.arange(azimuth_order) + .5) * 2 * np.pi / azimuth_order
    radius = np.sqrt(1 - mu ** 2)
    x = (radius[:, None] * np.cos(phi)[None, :]).ravel()
    y = (radius[:, None] * np.sin(phi)[None, :]).ravel()
    z = np.broadcast_to(-mu[:, None], (radial_order, azimuth_order)).ravel()
    limb = 1 - .45 * (1 - mu) - .2 * (1 - mu) ** 2
    weights = np.broadcast_to((weights / 2 * mu * limb)[:, None] * 2 * np.pi / azimuth_order, (radial_order, azimuth_order)).ravel()
    return x, y, z, weights


def flux_series(times, equatorial_period=10, inclination_degrees=60, latitude_degrees=20, shear=0, lifetime_days=np.inf, peak_day=50, longitudes=(0,), order=(64, 240)):
    x, y, z, quiet_weight = disk_quadrature(*order)
    inclination = np.deg2rad(inclination_degrees)
    latitude = np.deg2rad(latitude_degrees)
    local_period = equatorial_period / (1 - shear * np.sin(latitude) ** 2)
    result = []
    for time in times:
        contrast = .7 * np.exp(-.5 * ((time - peak_day) / lifetime_days) ** 2)
        brightness = np.ones(len(x))
        for initial_longitude in longitudes:
            phase = initial_longitude + 2 * np.pi * time / local_period
            direction = [np.cos(latitude) * np.sin(phase), np.sin(latitude) * np.sin(inclination) - np.cos(latitude) * np.cos(phase) * np.cos(inclination), -np.sin(latitude) * np.cos(inclination) - np.cos(latitude) * np.cos(phase) * np.sin(inclination)]
            region = np.exp(-(1 - x * direction[0] - y * direction[1] - z * direction[2]) / np.deg2rad(8) ** 2)
            brightness *= 1 - contrast * region
        result.append(np.sum(quiet_weight * brightness) / np.sum(quiet_weight))
    return np.array(result), local_period


def sinusoidal_fit_power(times, values, frequencies):
    """Floating-mean least-squares fraction of variance explained, not a FAP."""
    centered_time = np.asarray(times) - np.mean(times)
    values = np.asarray(values)
    baseline = np.sum((values - np.mean(values)) ** 2)
    result = []
    for start in range(0, len(frequencies), 256):
        phase = 2 * np.pi * np.asarray(frequencies[start:start + 256])[:, None] * centered_time[None, :]
        design = np.stack((np.ones_like(phase), np.sin(phase), np.cos(phase)), axis=-1)
        normal = np.einsum('fni,fnj->fij', design, design)
        right = np.einsum('fni,n->fi', design, values)
        fitted = np.linalg.solve(normal, right[..., None])[..., 0]
        residual = values[None, :] - np.einsum('fni,fi->fn', design, fitted)
        result.extend(1 - np.sum(residual ** 2, axis=1) / baseline)
    return np.array(result)