"""
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 sys
import numpy
from typing import Dict
from . import utils
import logging
import pkg_resources
logger = logging.getLogger(__name__)


def wd_stream(radial_angle: numpy.array, wndir: numpy.array) -> numpy.ndarray:
    wd = 6 * numpy.cos(radial_angle - wndir) * 10**(-2)
    return wd


class WaveDoppler():
    '''Class wave_doppler defines the computation of the wave doppler that is 
    seen by the instrument for STREAM and SKIM concepts
    Args:
         parameters: Simulation settings
    '''

    def __init__(self, parameters):
        # Seed is not used ? 
        self.nseed = parameters.nseed + 2
        self.type_satellite = parameters.type_satellite
        self.beam_angle = parameters.list_angle

    def generate(self, wind_norm: numpy.array, wind_dir: numpy.array,
                 windr: numpy.array, mss: numpy.array, hs: numpy.array,
                 ice: numpy.array, sigma0: numpy.array,
                 radial_angle: numpy.array, uuss: numpy.array,
                 vuss: numpy.array) -> Dict[str, numpy.array]:
        ''' Build errors corresponding to the wave Doppler seen by the
        instrument on the satellite swath.
        Args:
            radial_angle: Angle refered to the across track direction
            wind_dir: Wind direction refered to the East
            wind_norm: Wind velocity
            mss: Mean square slope (mssu + mssc or mssclose)
            hs: significant wave height
            ice: ice
            sigma0: sigma0
            uuss: stokes drift Eastward componant
            vuss: stokes drift Northward componant
            radial_angle: angle refered to the East
        '''
        if self.type_satellite == 'stream':
            uwd = wd_stream(radial_angle, wind_dir)
            ussr = numpy.full(numpy.shape(uwd), numpy.nan)
        elif self.type_satellite == 'skim':
            ussr = utils.proj_radial(uuss, vuss, radial_angle)
            uwd_all = []
            _comp = compute_wd_ai_par
            for b in range(len(self.beam_angle)):
                uwd = _comp(wind_norm[:, b], windr[:, b], mss[:, b],
                                hs[:, b], ice[:, b], sigma0[:, b],
                                ussr[:, b], radial_angle[:, b],
                                self.beam_angle[b])
                uwd_all.append(uwd)
            uwd_all = numpy.transpose(numpy.array(uwd_all))
        else:
            logger.error(f'Unknown concept {self.type_satellite}')
            sys.exit(1)
        return {"wd": uwd_all, "ur_uss": ussr}


    def correct(self, sgrid:dict, dic_input:dict, wind_norm: numpy.array,
                windr:numpy.array, ice: numpy.array,
                sigma0: numpy.array, radial_angle: numpy.array,
                dazim: float)-> Dict[str, numpy.ndarray]:
        if self.type_satellite == 'stream':
            return {}
        elif self.type_satellite == 'skim':
            from . import simulate_spectrum
            _sim = simulate_spectrum.simulate_spect_skim
            usr_est, mss_est, hs_est = _sim(sgrid, dic_input, self.beam_angle,
                                            radial_angle, delta_azim=dazim,
                                            fc=1, nseed=self.nseed)
            uwd_all = []

            _comp = compute_wd_ai_par
            for b in range(len(self.beam_angle)):
                uwd = _comp(wind_norm[:, b], windr[:, b], mss_est[:, b],
                            hs_est[:, 0], ice[:, b], sigma0[:, b],
                            usr_est[:, b], radial_angle[:, b],
                            self.beam_angle[b])
                uwd_all.append(uwd)
            uwd_all = numpy.transpose(numpy.array(uwd_all))
            usr_est = usr_est
            return {"wd_est": uwd_all, "ur_uss_est": usr_est}

def compute_wd_old_par(uuss: numpy.array, vuss:numpy.array, ucur: numpy.array,
                       vcur: numpy.array, uwnd: numpy.array, vwnd: numpy.array,
                       radial_angle:numpy.array, beam_angle: float
                       ) -> Dict[str, numpy.ndarray]:
    ''' Compute Wave doppler for SKIM instrument  using old parametrisation
    uwd = Gr * ussr and Gr = a * log(b + c/nwr) * (1- tanh(angle))
    Args:
        uuss: Stokes drift eastward componant
        vuss: Stokes drift northward componant
        ucur: Current eastward componant
        vcur: Current northward componant
        uwnd: Wind eastward componant
        vwnd: Wind northward componant
        radial_angle: Angle refered to the East
        beam_angle: View angle of beam
    '''
    ur_uss = utils.proj_radial(uuss, vuss, radial_angle)
    nwr = numpy.sqrt((uwnd - ucur)**2 + (vwnd - vcur)**2)
    nwr[nwr==0] = numpy.nan
    _angle = numpy.deg2rad((beam_angle - 25) / 10)
    GR = 25 * (0.82 * numpy.log(0.2 + 7/nwr)) * (1 - numpy.tanh(_angle))
    GP = 0
    uwd_noerr = GR * ur_uss
    cshape = numpy.shape(uwd_noerr)
    noise = numpy.random.normal(0, abs(uwd_noerr) * 0.25, cshape[0])
    wd = uwd_noerr + noise
    return {"ur_uss": ur_uss, "wd_noerr": uwd_noerr, "wd": wd}


def compute_wd_ai_par(nwnd: numpy.array, wndr: numpy.array, mss: numpy.array,
                      hs: numpy.array, ice: numpy.array, sigma0: numpy.array,
                      ussr: numpy.array, radial_angle: numpy.array,
                      beam_angle: float) -> Dict[str, numpy.array]:
    ''' Compute wave doppler using coefficients learned from ww3 data 
    Args:
        nwmd: Wind velocity
        wndr: Wind direction refered to the East
        mss: Mean square slope (mssu + mssc or mssclose)
        hs: significant wave height
        ice: ice
        sigma0: sigma0
        ussr: stokes radial componant
        radial_angle: angle refered to the East
        beam_angle: beam view incidence angle
    '''
    # Load Coefficents
    coeff_path = pkg_resources.resource_filename('skimulator',
                                                 'share/coeffr.npy')
    coeff = numpy.load(coeff_path, allow_pickle=True)[()]
    coeff_f1ur = coeff['f1UWDRglob{:d}deg'.format(int(beam_angle))]
    coeff_b1ur = coeff['b1UWDRglob{:d}deg'.format(int(beam_angle))]
    coeff_m = coeff['x_4dmean_{:d}deg'.format(int(beam_angle))]
    coeff_std = coeff['x_4dstd_{:d}deg'.format(int(beam_angle))]
    xlabel = coeff['xx{:d}'.format(int(beam_angle))]
    ncoeffur = 7 #len(coeff_4dm)
    # -- Compute norm --
    # Construct matrix of input data
    cshape = numpy.shape(wndr)
    ussr_max = Usr_95(wndr)
    ussr_min = Usr_5(wndr)
    diff_max = ussr_max - ussr
    diff_min = ussr_min - ussr
    ussr[numpy.where(diff_max < 0)] = ussr_max[numpy.where(diff_max < 0)]
    ussr[numpy.where(diff_min > 0)] = ussr_min[numpy.where(diff_min > 0)]
    wndr_min = + wndr
    wndr_min[numpy.where(abs(wndr) >= 8.)] = 8.
    nwnd_min = + nwnd
    nwnd_min[numpy.where(nwnd >= 8.)] = 8.
    mss_max = mss_95(nwnd)
    mss_min = mss_5(nwnd)
    diff_max = (mss_max - mss)
    diff_min = (mss_min - mss)
    mss[numpy.where(diff_max < 0)] = mss_max[numpy.where(diff_max < 0)]
    mss[numpy.where(diff_min > 0)] = mss_min[numpy.where(diff_min > 0)]

    #    noise_nwnd = numpy.random.normal(0, 10, cshape[0])
    #    wndr = wndr + noise_nwnd
    _ind = numpy.where(mss <= 8e-4)
    mss[_ind] = 8e-4

    mat_noerr = numpy.full((cshape[0], ncoeffur), numpy.nan)
    mat_noerr[:, 0] = ussr
    mat_noerr[:, 1] = wndr
    # mat_noerr[:, 2] = numpy.sign(wndr) *  wndr_min
    mat_noerr[:, 2] = nwnd
    mat_noerr[:, 3] = nwnd_min
    mat_noerr[:, 4] = hs
    mat_noerr[:, 5] = 1. / mss
    mat_noerr[:, 6] = 1. / (mss + numpy.log(nwnd + 0.7) * 0.009)

    # Compute cross product
    cross = utils.cross_product(mat_noerr, ncoeffur, cshape[0])
    ncoeffcross = len(coeff_m)
    # Normalize with coefficents
    for i in range(ncoeffcross):
        cross[:, i] = (cross[:, i] - coeff_m[i]) / coeff_std[i]

    shape_b1 = numpy.shape(coeff_b1ur)[0]
    Uwd = utils.reconstruct_var(ncoeffcross, cshape[0], coeff_f1ur,
                                    coeff_b1ur, cross, xlabel)
    if beam_angle == 6:
        sigma_ice = 2.5
    else:
        sigma_ice = 1.
    Uwd = Uwd * (1 - sigma_ice * ice / sigma0)
    Uwd[_ind] = 0.001
    Uwd[numpy.where(ice==1)] = 0.001
    return Uwd


def Usr_95(wndr):
    # wndr : radial wind speed [m/s]
    # provides 95th percentile of USR (max tolerable)
    wndr = wndr - 0.75 # middle point to lower bound
    mawndr = numpy.maximum(wndr, -8)
    miwndr = numpy.minimum(wndr, 8)
    res = 0.0200072*wndr - 0.0073556*miwndr - 0.00290607*mawndr + 0.03929183
    return res


def Usr_5(wndr):
    # wndr : radial wind speed [m/s]
    # provides 95th percentile of USR (min tolerable)
    wndr = wndr - 0.75 # middle point to lower bound
    mawndr = numpy.maximum(wndr, -8)
    miwndr = numpy.minimum(wndr, 8)
    res = 0.01896863*wndr - 0.00267442*miwndr - 0.00657048*mawndr - 0.02238212
    return res


def mss_5(wnd):
    # wnd : wind speed [m/s]
    # provides 5th percentile of mss (min tolerable)
    wnd = wnd - 0.5 # middle point to lower bound
    mawnd = numpy.maximum(wnd, 5)
    miwnd = numpy.minimum(wnd, 20)
    res = (0.00152612277*wnd - 4.73687795e-06*mawnd + 2.14390064e-04*wnd**2
           - 2.28370767e-04*wnd*mawnd)
    return res


def mss_95(wnd):
    # wnd : wind speed [m/s]
    # provides 95th percentile of mss (max tolerable)
    mawnd = numpy.maximum(wnd, 5)
    miwnd = numpy.minimum(wnd, 20)
    wnd = wnd - 0.5 # middle point to lower bound
    res = (2.71004327e-03 * wnd + 1.27761483e-03*mawnd - 1.68132538e-03*miwnd
           - 4.400380528e-5*wnd**2 - 1.35712627e-05)
    return res
