"""Fixed aperture/astrometric selections, never a photometric cleaning of binaries."""
import numpy as np


def angular_distance(ra, dec, centre=(132.85, 11.81)):
    a, d = np.radians(ra), np.radians(dec)
    a0, d0 = np.radians(centre)
    haversine = np.sin((d-d0)/2)**2 + np.cos(d)*np.cos(d0)*np.sin((a-a0)/2)**2
    return np.degrees(2*np.arcsin(np.sqrt(np.clip(haversine, 0, 1))))


def select_members(columns):
    quality = np.ones(len(columns['ra']), dtype=bool)
    for name in ['parallax', 'pmra', 'pmdec', 'phot_g_mean_mag', 'phot_bp_mean_mag', 'phot_rp_mean_mag']:
        quality &= np.isfinite(columns[name])
    for name in ['parallax_error', 'pmra_error', 'pmdec_error']:
        quality &= np.isfinite(columns[name]) & (columns[name] > 0)
    quality &= columns['visibility_periods_used'] >= 9
    quality &= (columns['phot_g_mean_flux_over_error'] > 50)
    quality &= (columns['phot_bp_mean_flux_over_error'] > 10)
    quality &= (columns['phot_rp_mean_flux_over_error'] > 10)
    # These are fixed literature-inspired cuts, not membership probabilities.
    pm_radius = np.hypot(columns['pmra']+11, columns['pmdec']+2.9)
    plx_difference = np.abs(columns['parallax']-1.15)
    selected = quality & (pm_radius < .7) & (plx_difference < .12)
    # A separately stated error-aware sensitivity selection includes correlations.
    errors = np.column_stack([columns[name] for name in ['parallax_error', 'pmra_error', 'pmdec_error']])
    residuals = np.column_stack([columns['parallax']-1.15, columns['pmra']+11, columns['pmdec']+2.9])
    pulls = np.full(len(quality), np.nan)
    for index in np.flatnonzero(quality):
        covariance = np.diag(errors[index]**2)
        for a, b, name in [(0, 1, 'parallax_pmra_corr'), (0, 2, 'parallax_pmdec_corr'), (1, 2, 'pmra_pmdec_corr')]:
            covariance[a, b] = covariance[b, a] = columns[name][index]*errors[index, a]*errors[index, b]
        covariance += np.diag([.05**2, .2**2, .2**2])
        if not np.all(np.linalg.eigvalsh(covariance) > 0):
            raise ValueError('Non-positive astrometric covariance')
        pulls[index] = residuals[index] @ np.linalg.solve(covariance, residuals[index])
    sensitivity = quality & (pulls < 16.27)
    radius = angular_distance(columns['ra'], columns['dec'])
    return {'quality': quality, 'fixed': selected, 'covarianceSensitivity': sensitivity,
            'radiusDegrees': radius, 'pmRadiusMasYr': pm_radius, 'parallaxDifferenceMas': plx_difference}