"""Finite brightness-basis visibility operators, with explicit numerical resolution."""
import numpy as np


def operator(phases, latitudes, longitudes, inclination, eclipse):
    latitude_sine = -1 + (np.arange(latitudes)+.5)*2/latitudes
    latitude_cosine = np.sqrt(1-latitude_sine**2)
    longitude = -np.pi + (np.arange(longitudes)+.5)*2*np.pi/longitudes
    stellar = latitude_cosine[:, None]*np.cos(longitude)[None, :]
    east = latitude_cosine[:, None]*np.sin(longitude)[None, :]
    north = np.broadcast_to(latitude_sine[:, None], stellar.shape)
    basis = np.array([np.ones(longitudes)] + [function(order*longitude) for order in range(1, 7) for function in [np.cos, np.sin]])
    matrix = []
    separation, stellar_radius = 4.882274180849331/.1577, 1/.1577
    observer_latitude = np.pi/2-inclination
    for phase in phases:
        observer_longitude = np.pi-2*np.pi*phase
        projection = (stellar*np.cos(observer_latitude)*np.cos(observer_longitude)
                      + east*np.cos(observer_latitude)*np.sin(observer_longitude) + north*np.sin(observer_latitude))
        weight = np.maximum(projection, 0.)
        if eclipse and np.cos(2*np.pi*phase) < 0:
            # Screen right is global +x; projected east points left at secondary.
            x = np.sin(observer_longitude)*stellar - np.cos(observer_longitude)*east
            y = (-np.sin(observer_latitude)*np.cos(observer_longitude)*stellar
                 - np.sin(observer_latitude)*np.sin(observer_longitude)*east + np.cos(observer_latitude)*north)
            star_x = separation*np.sin(2*np.pi*phase)
            star_y = separation*np.cos(inclination)*np.cos(2*np.pi*phase)
            weight = np.where((x-star_x)**2+(y-star_y)**2 < stellar_radius**2, 0., weight)
        longitude_weight = np.sum(weight, axis=0)*4*np.pi/(latitudes*longitudes)
        matrix.append(np.einsum('ml,l->m', basis, longitude_weight))
    return np.array(matrix)


def observed_window_rank(times, epoch, period, source_indices):
    phase = ((times-epoch)/period)%1
    # Preserve actual sample phases: every 32nd integration and all secondary contacts.
    selected = (source_indices%32 == 0) | (np.abs(phase-.5) < .045)
    selected &= np.minimum(phase, 1-phase) > .045
    phases = phase[selected]
    inclination = np.deg2rad(82.27012826519857)
    rows = []
    for latitudes, longitudes in [(24, 256), (48, 512)]:
        for label, angle, eclipse in [('edge_on_phase_only', np.pi/2, False), ('observed_inclination_phase_only', inclination, False), ('observed_inclination_and_eclipse', inclination, True)]:
            matrix = operator(phases, latitudes, longitudes, angle, eclipse)
            singular = np.linalg.svd(matrix, compute_uv=False)
            rows.append(dict(kind=label, latitude_bands=latitudes, longitude_cells=longitudes,
                             singular_values_normalized=(singular/singular[0]).tolist()))
    return dict(kind='Finite latitude-independent intensity basis through order six; not an arbitrary two-dimensional map',
                sample_count=int(selected.sum()), actual_source_indices=source_indices[selected].tolist(),
                selection='Actual timestamps: every 32nd original integration index plus all phase within 0.045 cycle of secondary; primary omitted.',
                results=rows,
                interpretation='Edge-on odd-mode zeros have finite quadrature residue. Compare both resolutions before interpreting small singular values. Eclipse can add sensitivity to phase-curve null modes; finite-basis rank is not uniqueness of a general map.')