"""
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/>.
"""
""" Error due to Sigma gradient within a beam simulation"""

import sys
import numpy
from typing import Dict, Tuple
from . import utils
import logging
logger = logging.getLogger(__name__)


def compute_pdf_dsigma(beam_angle: float)-> Tuple[float, float]:

    # Compute sigma for pdf altika
    inc = numpy.deg2rad(beam_angle)
    sig_altika = numpy.arange(7, 15)
    siglin = 10.**(sig_altika / 10)
    R = 0.55
    mss = R**2 / siglin
    dmssdsig0 = -R**2 / siglin**2

    xe = numpy.linspace(-0.02, 0.02, 1001)
    dxe = xe[1] - xe[0]

    # Compute sum of pdfs
    a = 19.5
    nu = 4.635 * numpy.exp(-0.03 * (sig_altika - 12)**2)
    nu_m = numpy.array([nu, ]* len(xe))
    c = 0.000045 / (10**((sig_altika + 1.5)**2 / 100)) + 4.9e-7 * sig_altika
    c_m = numpy.array([c, ]*len(xe))
    xe_m = numpy.array([xe, ] * len(sig_altika)).transpose()
    sigbeam = (R**2.*numpy.exp(-(numpy.tan(inc)**2 / mss))
                             / (mss * numpy.cos(inc)**4))
    dmss = 0.0001
    dsigdbeam = (R**2*numpy.exp(-(numpy.tan(inc)**2./(mss + dmss)))
                 / ((mss + dmss)*numpy.cos(inc)**4) - sigbeam) / dmss
    stretched_fac = siglin * dmssdsig0 * dsigdbeam / sigbeam
    stretched_fac_m =  numpy.array([stretched_fac, ]*len(xe))
    dxe = xe[1] - xe[0]

    mat_pdf = 10**((a - (3 + nu_m) * numpy.log(1 + xe_m**2
                    / ((1 + nu_m) * c_m)))/10)
    sum_pdf = numpy.array([numpy.sum(mat_pdf, axis=1), ] * len(mss)).transpose()
    norm_pdf0 = mat_pdf / (sum_pdf) # * dxe)
    stretched_pdf = abs(1/stretched_fac)*10**((a - (3 + nu_m) * numpy.log(1
                         + (xe_m/stretched_fac)**2
                         / ((1 + nu_m) * c_m)))/10) #/sum_pdf

    return xe, stretched_pdf


class DSigma():
    '''Class that simulates variation of sigma within a footprint. It is
    derived from Altika observation
    Args:
        parameters (nseed, sat_elev)

    '''
    def __init__(self, parameters):
        self.nseed = parameters.nseed
        self.sat_elev = parameters.sat_elev / 10**3
        self.beam_angle = parameters.list_angle
        self.vsat = parameters.vsat


    def generate(self, sigma0: numpy.ndarray,
                  ac_angle: numpy.ndarray)-> Dict[str, numpy.ndarray]:
        '''Generate statistical dsigma error within a footprint for one beam
        Args: 
            beam_angle: angle view from beam
            sigma0: sigma0 for one beam
            ac_angle: angle from across track direction
        '''
        dsigma_all = []
        for b in range(len(self.beam_angle)):
            xe, pdf = compute_pdf_dsigma(self.beam_angle[b])
            inc = numpy.deg2rad(self.beam_angle[b])
            # ratio altika - skim
            r = 0.66
            cshape = numpy.shape(sigma0[:, b])
            dsigma = numpy.full(cshape, numpy.nan)
            # to change, very slow
            for i in range(cshape[0]):
                if not numpy.isfinite(sigma0[i, b]):
                    continue
                if numpy.isnan(sigma0[i, b]):
                    continue

                siglog = 10 * numpy.log10(sigma0[i, b])
                ind_sig = numpy.floor(siglog) - 7
                ind_sig = max(ind_sig, 7)
                ind_sig = int(min(ind_sig, 0))
                pdf_stretched = pdf[:, ind_sig]
                distribution = pdf_stretched / numpy.sum(pdf_stretched)
                try:
                    res = numpy.random.choice(xe, size=1, p=distribution)
                    dsigma[i] = res[0]
                except:
                    dsigma[i] = numpy.nan
            dsigma = (r * dsigma / (10*numpy.log(10)) * self.sat_elev
                      * numpy.sin(inc))
            dsigma = utils.convert_dbkm2ms(dsigma, ac_angle[:, b],
                                           self.beam_angle[b], self.vsat)
            dsigma_all.append(dsigma)
        return {"dsigma": numpy.transpose(numpy.array(dsigma_all))}
