#double f_int(double t, int n)
#       {
#           return (pow(t,n))/(exp(t)-1)*(1-exp(-t))
#       }
#void _nlsfBlochGruenseimD(#double P1, double P2, double P3, double P4, double x,double& y)
#{
#       int n;
#       double y0;
#       n=2;
#       double t;
#       t=1e-6;
#       y0=36.6465;
#       y=y0;
#       double c;
#       double d;
#       double e;
#       double dIntegral;
#
#       do // loop over n
#       {
#           do //loop over t+dStep
#       {
#           // Initialization
#           double dIntegral = 0.0;
#           double dTrapezia = 0.0;
#           // Steps, or Precision.
#           double dStep = 1e-12;
#           // Trapezia area.
#           dTrapezia = 0.5 * (f_int(t, n)+ f_int((t+dStep),n)) * dStep;
#           // Accumulate area.
#           dIntegral += dTrapezia;     // intregral calculation.
#           t += dStep;
#           }while(t<=P1/x);
#       // Perform integrate by trapezoidal rule.
#
#           if(n=2)
#               c=P2;
#           else if(n=3)
#           d=P3;
#           else
#           e=P4;   //n=5
#           n=n+1;
#
#       }while (n<=5);
#
#           y = y + c * pow(x,2) * dIntegral+ d * pow(x,3) * dIntegral+e * pow(x,5) * dIntegral;
#   // End of editable part
#}
#
def f_int(t, n):
    return (t**n)/((math.exp(t)-1)*(1-math.exp(-t)))

def BD_Func(P1, P2, P3, P4, x, y):
    c, d, e =  None, None, None
    n = 2
    t = 1e-6
    y0 = 36.6465
    y = y0
    while n <= 5:
        while t <= (P1/x):
            dInt, dTrap = 0, 0
            dStep = 1e-5

            #Trap. Area
            dTrap = 0.5 * ((f_int(t, n)+ f_int((t+dStep),n)) * dStep)
            dInt += dTrap
            t += dStep
        if n == 2:
            c = P2
        elif n == 3:
            d = P3
        else:
            e = P4
        n += 1
    y = y + c * (x**2) * dIntegral + d * (x**3) * dIntegral+e * (x**5) * dIntegral
    return y