# -*- coding: utf-8 -*-
"""
Created on Thu Apr  2 15:00:25 2015

@author: adrienkuntz
"""

import numpy as np
cimport numpy as np
from scipy import integrate
from scipy import optimize
cimport cython_gsl as gsl
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))
    
 
cdef double 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
    
    
cdef double [::1] chitab = np.linspace(0., chiCMB, chinum)
cdef int i

cdef double [::1] Wktab = np.zeros(chinum)
cdef double [::1] Wgtab = np.zeros(chinum)
cdef double [::1] redshifttab = np.zeros(chinum)

for i in range(chinum) :
    redshifttab[i] = redshift(chitab[i])
       

for i in range(chinum) :
    Wktab[i] = Wk(chitab[i])
    Wgtab[i] = Wg(chitab[i])
    redshifttab[i] = redshift(chitab[i])
    
Wkfile = open('Wk.dat', 'w')
for i in range(chinum) :
    Wkfile.write('{} {}\n'.format(chitab[i], Wktab[i]))
Wkfile.close()

Wgfile = open('Wg.dat', 'w')
for i in range(chinum) :
    Wgfile.write('{} {}\n'.format(chitab[i], Wgtab[i]))
Wgfile.close()


cdef gsl.gsl_interp_accel *accWk = gsl.gsl_interp_accel_alloc()
cdef gsl.gsl_interp *interpWk = gsl.gsl_interp_alloc(gsl.gsl_interp_cspline, chinum) 
gsl.gsl_interp_init(interpWk, &chitab[0], &Wktab[0], chinum)

cdef gsl.gsl_interp_accel *accWg = gsl.gsl_interp_accel_alloc()
cdef gsl.gsl_interp *interpWg = gsl.gsl_interp_alloc(gsl.gsl_interp_cspline, chinum) 
gsl.gsl_interp_init(interpWg, &chitab[0], &Wgtab[0], chinum)

cdef gsl.gsl_interp_accel *accz = gsl.gsl_interp_accel_alloc()
cdef gsl.gsl_interp *interpz = gsl.gsl_interp_alloc(gsl.gsl_interp_cspline, chinum) 
gsl.gsl_interp_init(interpz, &chitab[0], &redshifttab[0], chinum)


i=0
monfichier = open('Wk.dat', 'r')
iterable = monfichier.readlines()
monfichier.close()
for ligne in iterable :
    part = ligne.split(' ')
    Wktab[i] = float(part[1])
    i = i + 1
i = 0
monfichier = open('Wg.dat', 'r')
iterable = monfichier.readlines()
monfichier.close()
for ligne in iterable :
    part = ligne.split(' ')
    Wgtab[i] = float(part[1])
    i = i + 1

cpdef double Wki(double chi) :
    return gsl.gsl_interp_eval(interpWk, &chitab[0], &Wktab[0], chi, accWk)
    
cpdef double Wgi(double chi) :
    return gsl.gsl_interp_eval(interpWg, &chitab[0], &Wgtab[0], chi, accWg)

cpdef double redshifti(double chi) :
    return gsl.gsl_interp_eval(interpz, &chitab[0], &redshifttab[0], chi, accz)