import datetime as dt
import numpy as np
import re
from math import cos, sin, pi, floor
from .constants import SECS_IN_MIN, SECS_IN_HR, EARTH_RADIUS
from .lib.coordinates import LocalCoord
from .gps_time import GPSTime

# Altitude of Ionospheric-pierce-point
IPP_ALT = 6821000

def get_alpha_beta(rcv_pos, el):
    geocentric_alt = np.linalg.norm(rcv_pos)
    alpha = np.pi/2 + el
    beta = np.arcsin(geocentric_alt*np.sin(alpha)/IPP_ALT)
    return alpha, beta

def get_slant_delay(rcv_pos, az, el, sat_pos, time, freq, vertical_delay):
    alpha, beta = get_alpha_beta(rcv_pos, el)
    slant_delay = vertical_delay * ((1 - ((EARTH_RADIUS * np.sin(beta)) /
                                          (EARTH_RADIUS + 3.5e5))**2)**(-0.5))
    return slant_delay

def closest_in_list(lst, val, num=2):
  """
    Returns two (`num` in general) closest values of `val` in list `lst`
    """
  idxs = sorted(lst, key=lambda x: abs(x - val))[:num]
  return sorted(list(lst).index(x) for x in idxs)


def get_header_line(headr, proprty):
    """
    :param headr: the header of the RINEX-file
    :param proprty: string-like property to search for (e.g. 'delta-utc')
    :return: the string of the ``headr`` containing ``property``
    """
    pattern = re.compile(proprty, re.IGNORECASE)
    for d in headr:
        if pattern.search(d):
            return d


def get_header_body(file_path):
    """
    Opens `file_path`, reads file and returns header and body
    separated with "END OF HEADER"
    :param file_path: path to RINEX-like file
    :return: header, body (arrays of lines)
    """
    with open(file_path) as fd:
        data = fd.readlines()
        for j, d in enumerate(data):
            if "END OF HEADER" in d:
                header_end = j
                break
    return data[:header_end], data[header_end + 1:]


def get_int_from_header(hdr, seq):
    """
    Returns the first int from the line that contains `seq` of lines `hdr`.
    In fact, _header_ here may not be header of RINEX/IONEX, just some set of lines.
    """
    return int(get_header_line(hdr, seq).split()[0])

def compute_grid_lats_lons(data):
  grid = np.array([], dtype='uint16')
  lats = np.array([])
  for j, line in enumerate(data[1:]):
    if "LAT" in line:
      lat, lon1, lon2, dlon, h = (float(line[x:x + 6]) for x in range(2, 32, 6))
      lats = np.append(lats, lat)
      row_length = (lon2 - lon1) / dlon + 1  # total number of values of longitudes
      next_lines_with_numbers = int(np.ceil(row_length / 16))
      elems_in_row = [
        min(16, int(row_length - i * 16)) for i in range(next_lines_with_numbers)
      ]
      row = np.array([], dtype='int16')
      for i, elem in enumerate(elems_in_row):
        row = np.append(
          row,
          np.array(
            [int(data[j + 2 + i][5 * x:5 * x + 5]) for x in range(elem)],
            dtype='int16',
          ),
        )
      if len(grid) > 0:
        grid = np.vstack((grid, row))
      else:
        grid = np.append(grid, row)
  lons = np.linspace(lon1, lon2, int(row_length))
  return (grid, lats, lons)


class IonexMap:
  def __init__(self, exp, data1, data2):
    self.exp = exp
    self.t1 = GPSTime.from_datetime(dt.datetime(*[int(d) for d in data1[0].split()[:6]]))
    self.t2 = GPSTime.from_datetime(dt.datetime(*[int(d) for d in data2[0].split()[:6]]))
    assert self.t2 - self.t1 == SECS_IN_HR
    assert len(data1) == len(data2)

    self.max_time_diff = SECS_IN_MIN*30
    self.epoch = self.t1 + self.max_time_diff

    self.grid_TEC1, self.lats, self.lons = compute_grid_lats_lons(data1)
    self.grid_TEC2, self.lats, self.lons = compute_grid_lats_lons(data2)

  def valid(self, time):
    return abs(time - self.epoch) <= self.max_time_diff

  @staticmethod
  def find_nearest(lst, val):
    return (np.abs(lst - val)).argmin()

  def get_TEC(self, pos, time):
    """
        Returns TEC in a position `pos`  of ionosphere
        :param pos: (lat, lon) [deg, deg]
        :return:
        """
    if pos[0] in self.lats and pos[1] in self.lons:
      lat = self.find_nearest(self.lats, pos[0])
      lon = self.find_nearest(self.lons, pos[1])
      E = self.grid_TEC1[lat][lon] + self.grid_TEC2[lat][lon]
      return E
    lat_idxs = closest_in_list(self.lats, pos[0])
    lon_idxs = closest_in_list(self.lons, pos[1])
    lat0, lat1 = self.lats[lat_idxs[0]], self.lats[lat_idxs[1]]
    lon0, lon1 = self.lons[lon_idxs[0]], self.lons[lon_idxs[1]]
    dlat = lat1 - lat0
    dlon = lon1 - lon0
    p = float(pos[0] - lat0) / dlat
    q = float(pos[1] - lon0) / dlon

    (E00, E10), (E01, E11) = self.grid_TEC1[lat_idxs[0]:lat_idxs[1] + 1, lon_idxs[0]:lon_idxs[1] + 1]
    TEC_1 = ((1 - p) * (1 - q) * E00 + p * (1 - q) * E01 + (1 - p) * q * E10 + p * q * E11)
    (E00, E10), (E01, E11) = self.grid_TEC2[lat_idxs[0]:lat_idxs[1] + 1, lon_idxs[0]:lon_idxs[1] + 1]
    TEC_2 = ((1 - p) * (1 - q) * E00 + p * (1 - q) * E01 + (1 - p) * q * E10 + p * q * E11)

    return (1 - (time - self.t1)/SECS_IN_HR)*TEC_1 + ((time - self.t1)/SECS_IN_HR)*TEC_2

  def get_delay(self, rcv_pos, az, el, sat_pos, time, freq):
    # To get a delay from a TEC map, we need to calculate
    # the ionospheric pierce point, geometry described here
    # https://en.wikipedia.org/wiki/Ionospheric_pierce_point
    alpha, beta = get_alpha_beta(rcv_pos, el)
    conv = LocalCoord.from_ecef(rcv_pos)
    gamma = np.pi - alpha - beta
    geocentric_alt = np.linalg.norm(rcv_pos)
    ipp_dist = geocentric_alt*np.sin(gamma)/np.sin(beta)
    ipp_ned = conv.ecef2ned(sat_pos)*(ipp_dist)/np.linalg.norm(sat_pos)
    ipp_geo = conv.ned2geodetic(ipp_ned)
    factor = 40.30E16 / (freq**2) * 10**(self.exp)
    vertical_delay = self.get_TEC(ipp_geo, time) * factor
    slant_delay = get_slant_delay(rcv_pos, az, el, sat_pos, time, freq, vertical_delay)
    return slant_delay

  @staticmethod
  def round_to_grid(number, base):
    return int(base * round(float(number) / base))


def parse_ionex(ionex_file):
  """
    :param ionex_file: path to the IONEX file
    :return: TEC interpolation function `f( (lat,lon), datetime )`
    """
  header, body = get_header_body(ionex_file)

  exponent = get_int_from_header(header, "EXPONENT")
  maps_count = get_int_from_header(header, "MAPS IN FILE")
  # =============
  # Separate maps
  # =============
  map_start_idx = []
  map_end_idx = []

  for j, line in enumerate(body):
    if "START OF TEC MAP" in line:
      map_start_idx += [j]
    elif "END OF TEC MAP" in line:
      map_end_idx += [j]
  if maps_count != len(map_start_idx):
    raise LookupError("Parsing error: the number of maps in the header "
                      "is not equal to the number of maps in the body.")
  if len(map_start_idx) != len(map_end_idx):
    raise IndexError("Starts end ends numbers are not equal.")
  map_dates = []
  for i in range(maps_count):
    date_components = body[map_start_idx[i] + 1].split()[:6]
    map_dates.append(dt.datetime(*[int(d) for d in date_components]))

  maps = []
  iono_map = iono_map_prev = None
  for m in range(maps_count):
    iono_map_prev = iono_map
    iono_map = body[map_start_idx[m] + 1:map_end_idx[m]]
    if iono_map and iono_map_prev:
      maps += [IonexMap(exponent, iono_map_prev, iono_map)]
  return maps


def klobuchar(pos, az, el, time, iono_coeffs):
  """
    Details are taken from [5]: IS-GPS-200H, Fig. 20-4
    Note: result is referred to the GPS L₁ frequency;
    if the user is operating on the GPS L₂ frequency, the correction term must
    be multiplied by γ = f₂²/f₁¹ = 0.6071850227694382
    :param pos: [lat, lon, alt] in radians and meters
    """

  tow = time.tow
  if pos[2] < -1E3 or el < 0:
    return 0.0
  if len(iono_coeffs) < 8:
    return None

  # earth centered angle (semi-circle)
  psi = 0.0137 / (el / pi + 0.11) - 0.022

  # subionospheric latitude/longitude (semi-circle)
  phi = pos[0] / pi + psi * cos(az)
  if phi > 0.416:
    phi = 0.416
  elif phi < -0.416:
    phi = -0.416
  lam = pos[1] / pi + psi * sin(az) / cos(phi * pi)

  # geomagnetic latitude (semi-circle) */
  phi += 0.064 * cos((lam - 1.617) * pi)

  # local time (s)
  tt = 43200.0 * lam + tow
  tt -= floor(tt / 86400.0) * 86400.0  # 0<=tt<86400

  # slant factor
  f = 1.0 + 16.0 * pow(0.53 - el / pi, 3.0)

  # ionospheric delay
  amp = iono_coeffs[0] + phi * (iono_coeffs[1] + phi *
                                (iono_coeffs[2] + phi * iono_coeffs[3]))
  per = iono_coeffs[4] + phi * (iono_coeffs[5] + phi *
                                (iono_coeffs[6] + phi * iono_coeffs[7]))
  if amp < 0.0:
    amp = 0.
  if per < 72000.0:
    per = 72000.0
  x = 2.0 * pi * (tt - 50400.0) / per

  mul = 5E-9
  if abs(x) < 1.57:
    mul = (5E-9 + amp * (1.0 + x * x * (-0.5 + x * x / 24.0)))
  return 2.99792458E8 * f * mul
