import os
import math
import warnings

import numpy as np

from astropy.time import Time
from astropy import constants as const
from astropy import units

from visnav.algo import tools
from visnav.algo.model import SystemModel, Camera, Asteroid
from visnav.iotools import objloader

from visnav.settings import *


class DidymosSystemModel(SystemModel):
    def __init__(self, target_primary=True, hi_res_shape_model=False, use_narrow_cam=True, res_mult=1.0):
        # gives some unnecessary warning about "dubious year" even when trying to ignore it
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            min_time = Time('2023-01-01 00:00:00', scale='utc', format='iso')

        common_kwargs_worst = {
            'sensor_size': (2048 * 0.0022, 1944 * 0.0022),                       # diagonal: 6.2mm
            'quantum_eff': 0.30,
            'px_saturation_e': 2200,                                             # snr_max = 20*log10(sqrt(sat_e)) dB
            'lambda_min': 350e-9, 'lambda_eff': 580e-9, 'lambda_max': 800e-9,
            'dark_noise_mu': 40, 'dark_noise_sd': 6.32, 'readout_noise_sd': 15,  # dark_noise_sd should be sqrt(dark_noise_mu)
            'emp_coef': 1,                                                       # dynamic range = 20*log10(sat_e/readout_noise))
            'exclusion_angle_x': 55,
            'exclusion_angle_y': 90,
        }
        common_kwargs_best = dict(common_kwargs_worst)
        common_kwargs_best.update({
            'quantum_eff': 0.4,
            'px_saturation_e': 3500,
            'dark_noise_mu': 25, 'dark_noise_sd': 5, 'readout_noise_sd': 5,
        })
        common_kwargs = common_kwargs_best

        narrow_cam = Camera(
            2048 * res_mult,   # width in pixels
            1944 * res_mult,   # height in pixels
            7.7,    # x fov in degrees  (could be 6 & 5.695, 5.15 & 4.89, 7.7 & 7.309)
            7.309,  # y fov in degrees
            f_stop=5,       # TODO: put better value here
            point_spread_fn=0.50,   # ratio of brightness in center pixel
            scattering_coef=2e-10,  # affects strength of haze/veil when sun shines on the lens
            **common_kwargs
        )
        wide_cam = Camera(
            2048 * res_mult,   # width in pixels
            1944 * res_mult,   # height in pixels
            61.5,     # x fov in degrees
            58.38,  # y fov in degrees
            f_stop=5,   # TODO: put better value here
            point_spread_fn=0.35,   # ratio of brightness in center pixel
            scattering_coef=2e-10,  # affects strength of haze/veil when sun shines on the lens
            **common_kwargs
        )

        if target_primary:
            if use_narrow_cam:
                mission_id = 'didy1n'
                limits = (
                    3.5,    # min_distance in km
                    6.5,    # min_med_distance in km
                    10.5,   # max_med_distance in km
                    10.5,   # max_distance in km
                    45,     # min_elong in deg
                    min_time,  # min time instant
                )
            else:
                mission_id = 'didy1w'
                limits = (
                    1.0,  # min_distance in km
                    1.1,  # min_med_distance in km
                    6.0,  # max_med_distance in km
                    10.5,  # max_distance in km
                    45,  # min_elong in deg
                    min_time,  # min time instant
                )
        else:
            if use_narrow_cam:
                mission_id = 'didy2n'
                limits = (
                    1.0,    # min_distance in km
                    1.65,    # min_med_distance in km
                    5.3,   # max_med_distance in km
                    7.0,   # max_distance in km
                    45,     # min_elong in deg
                    min_time,  # min time instant
                )
            else:
                mission_id = 'didy2w'
                limits = (
                    0.135,  # min_distance in km
                    0.28,  # min_med_distance in km
                    1.3,  # max_med_distance in km
                    1.3,  # max_distance in km
                    45,  # min_elong in deg
                    min_time,  # min time instant
                )

        super(DidymosSystemModel, self).__init__(
            asteroid=DidymosPrimary(hi_res_shape_model=hi_res_shape_model) if target_primary
                else DidymosSecondary(hi_res_shape_model=hi_res_shape_model),
            camera=narrow_cam if use_narrow_cam else wide_cam,
            limits=limits,
        )
        self.mission_id = mission_id
        self.sc_model_file = os.path.join(DATA_DIR, 'apex-x1-2019-05-28.obj')


class DidymosPrimary(Asteroid):
    # from ast frame (axis: +y, up: +z) to spacecraft coords (axis: +x, up: +z)
    ast2sc_q = np.quaternion(1, 1, 0, 0).normalized()

    # based on itokawa, Tatsumi et al 2018
    # "Regolith Properties on the S-Type Asteroid Itokawa Estimated from Photometrical Measurements"
    HAPKE_PARAMS = [
        # J, brightness scaling
        15,

        # th_p, average surface slope (deg), effective roughness, theta hat sub p
        40,

        # w, single scattering albedo (w, omega, SSA), range 0-1
        0.57,

        # b, SPPF asymmetry parameter (sometimes g?),
        #   if single-term HG, range is [-1, 1], from backscattering to isotropic to forward scattering
        #   if two-term HG (c>0), range is [0, 1], from isotropic to scattering in a single direction
        0.35,

        # c, second HG term for a more complex SPPF. Range [0, 1], from forward scattering to backward scattering.
        0.56,

        # B_SH0, or B0, amplitude of shadow-hiding opposition effect (shoe). If zero, dont use.
        0.98,

        # hs, or h or k, angular half width of shoe (rad)
        math.radians(0.05),

        # B_CB0, amplitude of coherent backscatter opposition effect (cboe). If zero, dont use.
        0,

        # hc, angular half width of cboe
        0.005,

        # extra mode selection, first bit: use K or not
        # NOTE: K generally not in use as phase angle changes so little inside one image and exposure is adjusted to
        #       increase overall brightness
        1,
    ]

    # fitted to above hapke params using approx-lunar-lambert.py: match_ll_with_hapke
    LUNAR_LAMBERT_PARAMS = [
        1.0,
        -3.2261e-02, 8.4991e-05, 1.4809e-06, -7.7885e-09, 8.7950e-12, 7.7132e-01,
        0, 0, 0
    ]

    def __init__(self, hi_res_shape_model=False):
        super(DidymosPrimary, self).__init__()
        self.name = 'Didymos Primary'

        self.image_db_path = None

        # use ryugu model for this, ryugu ~162m diameter, ryugu-big ~772m diameter (Didy2 & Didy1)
        self.target_model_file = os.path.join(DATA_DIR, 'ryugu+tex-d1-16k.obj')
        self.hires_target_model_file = os.path.join(DATA_DIR, 'ryugu+tex-d1-400k.obj')
        self.hires_target_model_file_textures = True

        self.constant_noise_shape_model = {
            '' : os.path.join(DATA_DIR, 'ryugu+tex-d1-16k.nsm'),  # same as target_model_file but includes error estimate
            'lo': os.path.join(DATA_DIR, 'ryugu+tex-d1-4k.nsm'),  # 1/4 the vertices
            'hi': os.path.join(DATA_DIR, 'ryugu+tex-d1-1k.nsm'),  # 1/16 the vertices
        }

        self.sample_image_file = None
        self.sample_image_meta_file = None

        self.real_shape_model = objloader.ShapeModel(
            fname=self.hires_target_model_file if hi_res_shape_model else self.target_model_file)
        self.render_smooth_faces = False if hi_res_shape_model else True

        self.reflmod_params = {
            1: DidymosPrimary.LUNAR_LAMBERT_PARAMS, # REFLMOD_LUNAR_LAMBERT
            2: DidymosPrimary.HAPKE_PARAMS, # REFLMOD_HAPKE
        }

        # for cross section, assume spherical object
        self.max_radius = 470      # in meters, maximum extent of object from asteroid frame coordinate origin
        self.mean_radius = 775/2   # in meters
        self.mean_cross_section = math.pi * self.mean_radius ** 2

        # epoch for orbital elements, 2019-Apr-27.0 TDB
        self.oe_epoch = Time(2458600.5, format='jd')

        # orbital elements (from https://ssd.jpl.nasa.gov/sbdb.cgi?sstr=2065803#content)
        # reference: JPL 134 (heliocentric ecliptic J2000)
        self.eccentricity = .3840204901532592
        self.semimajor_axis = 1.644267950023408 * const.au
        self.inclination = math.radians(3.408560852149408)
        self.longitude_of_ascending_node = math.radians(73.20707998527304)
        self.argument_of_periapsis = math.radians(319.3188822767833)
        self.mean_anomaly = math.radians(124.6176776030496)

        # other, not used
        self.aphelion = 2.275700534134692 * const.au
        self.perihelion = 1.012835365912124 * const.au
        self.orbital_period = 770.1180709731267 * 24 * 3600  # seconds
        # self.true_anomaly = math.radians(145.5260853202137 ??)

        # rotation period 11.92h from https://ssd.jpl.nasa.gov/sbdb.cgi?sstr=2065803#content
        self.rot_epoch = Time('J2000')
        self.rotation_velocity = 2 * math.pi / (2.2593 * 3600)  # rad/s

        # asteroid rotation axis in equatorial coordinates
        ra, de, tpm = 310, -84, 0   # ra=lon, de=lat

        self.rotation_pm = math.radians(tpm)
        self.axis_latitude, self.axis_longitude = \
            tuple(map(math.radians, (de, ra) if USE_ICRS else \
                tools.equatorial_to_ecliptic(de * units.deg, ra * units.deg)))

        self.precession_cone_radius = None
        self.precession_period = None
        self.precession_pm = None

        self.set_defaults()


class DidymosSecondary(Asteroid):
    ast2sc_q = DidymosPrimary.ast2sc_q

    def __init__(self, hi_res_shape_model=False):
        super(DidymosSecondary, self).__init__()
        self.name = 'Didymos Secondary'

        self.image_db_path = None

        # use ryugu model for this, ryugu ~162m diameter, ryugu-big ~772m diameter (Didy2 & Didy1)
        self.target_model_file = os.path.join(DATA_DIR, 'ryugu+tex-d2-16k.obj')
        self.hires_target_model_file = os.path.join(DATA_DIR, 'ryugu+tex-d2-400k.obj')

        self.constant_noise_shape_model = {
            '' : os.path.join(DATA_DIR, 'ryugu+tex-d2-16k.nsm'),  # same as target_model_file but includes error estimate
            'lo': os.path.join(DATA_DIR, 'ryugu+tex-d2-4k.nsm'),  # 1/4 the vertices
            'hi': os.path.join(DATA_DIR, 'ryugu+tex-d2-1k.nsm'),  # 1/16 the vertices
        }

        self.sample_image_file = None
        self.sample_image_meta_file = None

        self.real_shape_model = objloader.ShapeModel(
            fname=self.hires_target_model_file if hi_res_shape_model else self.target_model_file)
        self.render_smooth_faces = False if hi_res_shape_model else True

        self.reflmod_params = {
            1: DidymosPrimary.LUNAR_LAMBERT_PARAMS, # REFLMOD_LUNAR_LAMBERT
            2: DidymosPrimary.HAPKE_PARAMS, # REFLMOD_HAPKE
        }

        # for cross section, assume spherical object
        self.max_radius = 105     # in meters, maximum extent of object from asteroid frame coordinate origin
        self.mean_radius = 163/2  # in meters, dims = [206, 158, 132]
        self.mean_cross_section = math.pi * self.mean_radius ** 2

        # Distance = 1180 (1160-1220) m
        # L1, L2 = 999.3, 1354.4 m

        # epoch for orbital elements, 2019-Apr-27.0 TDB
        self.oe_epoch = Time(2458600.5, format='jd')

        # orbital elements (from https://ssd.jpl.nasa.gov/sbdb.cgi?sstr=2065803#content)
        # reference: JPL 134 (heliocentric ecliptic J2000)
        self.eccentricity = .3840204901532592
        self.semimajor_axis = 1.644267950023408 * const.au
        self.inclination = math.radians(3.408560852149408)
        self.longitude_of_ascending_node = math.radians(73.20707998527304)
        self.argument_of_periapsis = math.radians(319.3188822767833)
        self.mean_anomaly = math.radians(124.6176776030496)

        # other, not used
        self.aphelion = 2.275700534134692 * const.au
        self.perihelion = 1.012835365912124 * const.au
        self.orbital_period = 770.1180709731267 * 24 * 3600  # seconds
        # self.true_anomaly = math.radians(145.5260853202137 ??)

        # rotation period 11.92h from https://ssd.jpl.nasa.gov/sbdb.cgi?sstr=2065803#content
        self.rot_epoch = Time('J2000')
        self.rotation_velocity = 2 * math.pi / (11.92 * 3600)  # rad/s

        # asteroid rotation axis in equatorial coordinates
        ra, de, tpm = 310, -84, 0   # ra=lon, de=lat

        self.rotation_pm = math.radians(tpm)
        self.axis_latitude, self.axis_longitude = \
            tuple(map(math.radians, (de, ra) if USE_ICRS else \
                tools.equatorial_to_ecliptic(de * units.deg, ra * units.deg)))

        self.precession_cone_radius = None
        self.precession_period = None
        self.precession_pm = None

        self.set_defaults()
