import numpy as np
from scipy.special import erf
from scipy import integrate
import matplotlib.pyplot as plt


def VarLgauss(lc, blks, d):

    scl = (blks / lc) ** 2
    return (scl ** -d) * (
        (np.sqrt(2 * np.pi * scl) * erf(np.sqrt(scl / 2)) + 2 * np.exp(-0.5 * scl) - 2)
        ** d
    )


def VarLgaussSimp(lc, blks, d):

    A = lc / blks  # lc/L
    B = np.sqrt(2 * np.pi)  # square root of 2*pi
    C = erf((blks / lc) / np.sqrt(8))  # erf( (L/lc) / square root of 2)

    return (A * B * C) ** d


def arg_exp(t, lc, blks, d):

    scl = (blks / lc) ** 2

    aux1 = np.pi * erf(np.sqrt(scl / (4 * t)))
    aux2 = np.sqrt(np.pi) * (1 - np.exp(-scl / (4 * t))) * np.sqrt(4 * t / scl)
    return t * np.exp(-t) * ((aux1 - aux2) ** d)


def VarLexp3d(lc, blks):  # ic=5.378669493723924333  para lc 16
    d = 3
    # a=1/64/np.pi
    t = np.arange(0.000000001, 50, 0.001)
    var = np.empty(0)
    for blk in blks:
        y = arg_exp(t, lc, blk, d)
        var = np.append(var, 64 * np.pi * ((lc / (2 * np.pi * blk)) ** d) * np.trapz(y))

    return var


def argVarLexp2d(lc, blk):
    scl = float(blk / (2 * lc))
    f = lambda y, x: np.exp(-1 * np.sqrt(x ** 2 + y ** 2))

    res = integrate.dblquad(
        f, -scl, scl, lambda x: -scl, lambda x: scl, epsabs=1.49e-8, epsrel=1.49e-8
    )  # 0,1,lambda x: 0, lambda x: 1)

    return ((lc / blk) ** 2) * res[0]


def VarLexp2d(lc, blks):
    # if lc==1.33:
    # 	blks=np.append(np.arange(1,2,0.1),blks[1:])
    res = np.empty(0)
    for blk in blks:
        res = np.append(res, argVarLexp2d(lc, blk))

    return res