# -*- coding: utf-8 -*-
"""
Created on Fri Mar 13 13:36:51 2015

@author: adrienkuntz
"""

import numpy as np
from scipy import integrate
from scipy import optimize
from Params import *
#from scipy.special import sici
#from scipy.optimize import fsolve
#import matplotlib.pyplot as plt
#import math



### comoving distance as a function of redshift

def E(z): return np.sqrt(omegam0 * (1 + z)**3 + omegalambda0)
    
def chi(z) : return cl/H0 * integrate.quad(lambda x : 1/E(x), 0., z, epsrel = 1e-6, epsabs = 0)[0]
    
def chiv(z) : return np.vectorize(chi)(z)
    
def redshift(chi_): return optimize.brentq(lambda x : chiv(x) - chi_, 0., zCMB)

def redshiftv(chi) : return np.array(map(redshift, chi))
    
 
chiCMB = chi(zCMB)


#define the kernels ; for convenience I take an analytical dN/dz

def Wk(chi) : return 3./2 * H0**2 / cl**2 * omegam0 * chi * (chiCMB - chi) * (1 + redshift(chi)) / chiCMB

def dNdz(z) : return alpha / z0 * np.exp(lambda_ * np.log(z/z0)) * np.exp(-np.exp(beta * np.log(z/z0)))
    
intgl = integrate.quad(dNdz, 0., zCMB, epsabs = 0, epsrel = 1e-6)[0]

def Wg(chi_) :
    za = redshift(chi_)
    f = chi_ * integrate.quad(lambda z : dNdz(z) / intgl * (chi(z)-chi_) / chi(z), za, zCMB, epsabs=0, epsrel=1e-3)[0]
    return dNdz(za) / intgl * H0 / cl * E(za) + 3/2. * omegam0 * (H0/cl)**2 * (1+za) * (5*sn-2) * f