#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Nov 27 16:32:02 2017

@author: jared's account
"""

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import scipy

X=np.linspace(0.001,1,1001)
X2=np.linspace(0.0001,15,10001)
def fx(x):
    return x**(-1./2)/(np.exp(x)+1)

plt.plot(fx(X)/X**(-1./2.))
plt.show()

def px(x):
    return x**(-1./2.)/2.

def bb(x):
    return x**3*np.exp(-x)/(1-np.exp(-x))

a=.1

def pbb(x):
    return a*np.exp(-a*x)
 

n=1000000
ran=np.random.rand(n)
x=ran**2
I1=np.sum(fx(x)/px(x))/n
print(I1)

I2=scipy.integrate.simps(fx(X),X,even='avg')
print(I2)

plt.plot(X2,bb(X2))
plt.show()
I3=scipy.integrate.simps(bb(X2),X2)
print(I3*15/np.pi**4)
xn=-(1./a)*np.log(1-ran)
I4=np.sum(bb(xn)/pbb(xn))/n
print(I4*15/np.pi**4)