"""Functions to simulate scans and maps."""


import numpy as np
import numpy.random as ra
import os
from astropy.io import fits
from astropy.table import Table, vstack
import astropy.units as u
from collections.abc import Iterable

from .io import mkdir_p, locations
from .utils import tqdm, jit
from astropy.coordinates import SkyCoord
from astropy.time import Time
from astropy import log

try:
    import matplotlib.pyplot as plt

    HAS_MPL = True
except ImportError:
    HAS_MPL = False

__all__ = ["simulate_scan", "save_scan", "simulate_map"]


DEFAULT_PEAK_COUNTS = 100
COUNTS_TO_K = 0.03
DEFAULT_CAL_TEMP = 5
DEFAULT_CAL_OFFSET = DEFAULT_CAL_TEMP / COUNTS_TO_K


summary_header = """
SIMPLE  =                    T / file does conform to FITS standard
BITPIX  =                    8 / number of bits per data pixel
NAXIS   =                    0 / number of data axes
EXTEND  =                    T / FITS dataset may contain extensions
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
HIERARCH BackendName = 'NULL    ' / Backend name
CREATOR = 'NULL    '           / Software (incl. version)
DATE-OBS= '2016-10-03T14:59:08.753' / Observation time
EQUINOX =                   0. / Equinox of RA, Dec
EXPTIME =                   0. / Total integration time (seconds)
FITSVER = 'V.1.11  '           / FITS version
LST     =                    0 / Local sidereal time
HIERARCH LogFileName = 'NULL    ' / Name of the log file
HIERARCH NUSEBANDS =         0 / Number of sections
OBJECT  = 'W51     '           / Target source name
OBSID   = 'NULL    '           / Observer or operator initials
PROJID  = 'NULL    '           / ProjectID
HIERARCH RESTFREQ1 =  22235.08 / Rest frequency (MHz)
HIERARCH RESTFREQ2 =  22235.08 / Rest frequency (MHz)
HIERARCH RESTFREQ3 =  22235.08 / Rest frequency (MHz)
HIERARCH RESTFREQ4 =  22235.08 / Rest frequency (MHz)
HIERARCH ReceiverCode = 'CCB     ' / Receiver name
HIERARCH RightAscension = 5.07757730974885 / Target right ascension (radians)
HIERARCH Declination = 0.253290907695677 / Target declination (radians)
SCANGEOM= 'NULL    '           / Scan geometry
SCANMODE= 'NULL    '           / Mapping mode
SCANTYPE= 'NULL    '           / Scan astronomical type
SCANXVEL=                   0. / Tracking rate (optional,OTF)
SWTCHMOD= 'NULL    '           / Switch mode
HIERARCH ScheduleName = 'NULL    ' / Name of the schedule
TELESCOP= 'SRT     '           / Telescope name
VDEF    = 'OP      '           / Radial velocity definition
VFRAME  = 'LSRK    '           / Radial velocity reference frame
VRAD    =                    0 / Radial velocity
WOBUSED =                    0 / Wobbler used?

"""

_REFERENCE_MJD = 57000.5


def _apply_spectrum_to_data(spec_func, counts, nbin, bw=1000):
    """

    Examples
    --------
    >>> res = _apply_spectrum_to_data(lambda x: np.ones(x.size), 4, 3)
    >>> np.allclose(res, [4., 4., 4.])
    True
    >>> res = _apply_spectrum_to_data(lambda x: np.ones(x.size), [4, 2], 3)
    >>> np.allclose(res, [[4., 4., 4.], [2, 2, 2]])
    True
    >>> _apply_spectrum_to_data(lambda x: np.ones(x.size), 4, 1)
    4
    """
    if nbin == 1:
        return counts
    single = False
    if not isinstance(counts, Iterable):
        counts = [counts]
        single = True
    counts = np.asarray(counts)
    df = bw / nbin
    freqs = np.arange(0, bw, df)
    single_spec = spec_func(freqs)
    spec = np.zeros((len(counts), len(freqs)))
    for i, c in enumerate(counts):
        spec[i, :] += c * single_spec

    if single:
        return spec[0]
    return spec


def _standard_source_spectrum(counts, nbin, bw=1000, sigma=1):
    def spec_func(f):
        f = f - bw / 2
        return np.exp(-(f ** 2) / (2 * sigma ** 2))

    return _apply_spectrum_to_data(spec_func, counts, nbin, bw)


def _standard_bkg_spectrum(counts, nbin, bw=1000):
    def spec_func(f):
        sp = 1 + 0.1 * np.sin(2 * np.pi * 5 / bw * f) * (1 - f / bw)
        sp -= 0.5 * f / bw
        return sp

    return _apply_spectrum_to_data(spec_func, counts, nbin, bw)


def create_summary(filename, key_dict=None):
    if key_dict is None:
        key_dict = {}
    header = fits.Header.fromstring(summary_header, sep="\n")
    for key, value in key_dict.items():
        header[key] = value

    primary_hdu = fits.PrimaryHDU(header=header)
    hdul = fits.HDUList([primary_hdu])
    hdul.writeto(filename, overwrite=True)
    return filename


def _is_number(x):
    """ "Test if a string or other is a number

    Examples
    --------
    >>> _is_number('3')
    True
    >>> _is_number(3.)
    True
    >>> _is_number('a')
    False
    """
    try:
        float(x)
        return True
    except (ValueError, TypeError):
        return False


def _default_flat_shape(x):
    """A flat shape.

    Examples
    --------
    >>> _default_flat_shape(4314)