"""
Copyright (C) 2017-2021 OceanDataLab
This file is part of skimulator.

skimulator is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

skimulator is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with skimulator.  If not, see <http://www.gnu.org/licenses/>.
"""
import numpy
import logging
from typing import Dict

logger = logging.getLogger(__name__)


def instr_stream(radial_angle: numpy.ndarray, ac_angle: numpy.ndarray,
                 wind_norm: numpy.ndarray, wind_dir:numpy.ndarray,
                 cycle: float)-> numpy.ndarray:

    # - Compute NRCS from wind:
    diff_angle = radial_angle - wind_dir
    nrcs = (-47.5 + 17.5 * (wind_norm - 1)**0.3
            - numpy.sqrt(2)/2 *(1 - abs(numpy.cos(2 * diff_angle)))
            * (4 + 0.1* wind_norm))

    # NRCS at 0db is -29 db:
    snr = nrcs + 29
    instr = numpy.full(numpy.shape(wind_norm), numpy.nan)

    # Decifer formulae depending on SNR
    snr2 = + snr
    snr[snr==0] = 1
    N = 24 / (snr * 3.7 + 12)
    al_angle = ac_angle + numpy.pi/2
    instr = (abs(numpy.sin(al_angle)) * (10 + 4*numpy.tanh(snr/6))
             + (4 + 12 * (1 - numpy.tanh(snr / 6)))
             * abs(numpy.sin(al_angle))**N
             + 17 * (1 - numpy.tanh(snr / 6)))* 10**(-2)
    instr[snr>19] = (14.5 * abs(numpy.sin(al_angle[snr>19]))
                     + 4 * abs(numpy.sin(al_angle[snr>19]))**0.2
                     + 0.4) * 10**(-2)
    instr[snr2==0] = (abs(numpy.sin(al_angle[snr2==0])*10)
                      + 16 * abs(numpy.sin(al_angle[snr2==0]))**2
                      + 17) * 10**(-2)
    # - Knowledge of HH on top of VV will decrease error (about 20%)
    # Formulaes have been computed for 64 pulses, rms_instr factor will adapt
    # noise if number of pulses changes (sqrt(factor_of_pulses))
    instr = 0.8 * instr * numpy.sqrt(64/12000/2/cycle)
    return instr


def instr_skim_old(file_rms_instr: str, radial_angle: numpy.ndarray,
               rms_instr_factor:float):
    # - Load instrumental rms file into dictionnary
    theta, rms = numpy.loadtxt(file_rms_instr, usecols=(0, 1), unpack=True)
    dic_rms = {}
    with open(file_rms_instr) as f:
        for line in f:
            (key, val) = line.split()
            dic_rms[int(key)] = float(val)

    # - Compute rms corresponding to each radial angle
    # Refer radial angle to along track direction
    radial_angle_along = numpy.degrees((radial_angle + math.pi / 2)
                                       % math.pi)
    # Interpolate radial_angle with theta values
    step = theta[1]-theta[0]
    if((theta[1:]-theta[:-1]) != step).any():
        logger.error('instrumental rms file noise has not a constant step')
        sys.exit(1)
    thetaprev = numpy.floor(radial_angle_along/step)
    fac = radial_angle_along - thetaprev
    rms_theta = radial_angle * 0
    i = 0
    for itheta, ifac in zip(thetaprev, fac):
        keytheta = list(dic_rms.keys())[int(itheta)]
        keythetanext = list(dic_rms.keys())[int(itheta + step) % 180]
        rms_theta[i] = (dic_rms[keytheta] * ifac / step
                        + (step - ifac) / step * dic_rms[keythetanext])
        i += 1
    # - Errors array are the same size as the swath size
    instr = rms_theta * rms_instr_factor * 10**(-2)
    return instr


def instr_skim(beam_angle:float, ac_angle: numpy.ndarray, sigma0: numpy.array,
               cycle: float, instr_configuration: str)-> numpy.ndarray:
    if beam_angle == 12:
        if instr_configuration == 'A':
            _co = (-15.998, 0.657, -4.174, 15.260)
        if instr_configuration == 'B':
            _co = (16.409, 0.963, -6.232, 2.823)
        coeff = _co[0] * numpy.sin(ac_angle * _co[1] + _co[2])**2 + _co[3]
        sigma_ref = 10**(0.66)
    elif beam_angle == 6:
        if instr_configuration == 'A':
            _co = (17.272, -1.105, -2.977, 6.108)
        if instr_configuration == 'B':
            _co = (12.299, 1.667, -1.045, 14.071)
        coeff = _co[0] * numpy.sin(ac_angle * _co[1] + _co[2]) + _co[3]
        #igma_ref = 9.9
        sigma_ref = 10**(0.99)
    else:
        logger.error('Unknown instrumental parametrisation for {}'
                     ' angle'.format(beam_angle))
    coeff_random = (coeff * 10**(-2) * sigma0 / sigma_ref)
    if cycle < 0.0368:
        coeff_random = coeff_random * numpy.sqrt(0.0368/cycle)
    coeff_random[sigma0 < 10**-5] = numpy.nan
    return coeff_random


class Instrument():
    '''Class instrument defines the instrumental error for STREAM or SKIM
    concepts
    Args:
        parameters: Simulation settings (nseed, type_satellite,
     instr_configuration and cycle parameters)
    '''
    def __init__(self, parameters):
        self.nseed = parameters.nseed + 1
        self.type_satellite = parameters.type_satellite
        self.cycle = parameters.cycle
        self.instr_configuration = parameters.instr_configuration
        self.beam_angle = parameters.list_angle


    def generate(self, radial_angle: numpy.array, angle: numpy.array,
                 wind_norm: numpy.array, wind_dir: numpy.array,
                 sigma0: numpy.array) -> Dict[str, numpy.ndarray]:
        ''' Build errors corresponding to the instrument on the satellite
        swath.
        Args:
            radial_angle: Angle refered to the East
            angle: Angle refered to the across track direction
            wind_norm: wind velocity
            wind_dir: wind direction refered to the East
        '''
        instr_all = []
        for b in range(len(self.beam_angle)):
            if self.type_satellite == 'stream':
                instr = instr_stream(radial_angle[:, b], angle[:, b],
                                     wind_norm[:, b], wind_dir[:, b],
                                     self.cycle)
            elif self.type_satellite == 'skim':
                instr = instr_skim(self.beam_angle[b], angle[:, b],
                                   sigma0[:, b],
                                   self.cycle, self.instr_configuration)
            else:
                logger.error(f'Unknown concept {self.type_satellite}')
                sys.exit(1)
        # Seed to make reproducible noise
            numpy.random.seed(self.nseed)
            error = numpy.random.normal(0.0, abs(instr), numpy.shape(instr)[0])
            instr_all.append(error)
        return{"instr": numpy.transpose(numpy.array(instr_all))}
