''' This file contains the routine to calculate
	J_{\alpha, \beta, l}(k), as appears in 2.21 of the paper.

	It is the orginal FAST-PT code and has now been replaced by
	FASTPT.py.

	J. E. McEwen (c) 2016
'''
from __future__ import division
import numpy as np
from numpy.fft import fft, ifft , rfft, irfft , fftfreq
from numpy import exp, log, log10, cos, sin, pi, cosh, sinh , sqrt
from scipy.special import gamma
import sys
from time import time
from .fastpt_extr import p_window, c_window, pad_left, pad_right
from .gamma_funcs import g_m_vals, gamsn


log2=log(2.)


def check_conjugate(A):
	'''this function was only used for debugging.
	It is not important for actually running FAST-PT
	This module is used to check that signal is real,
	i.e. A(-omega)=A(omega)^*
	must be odd dimensional input
	'''

	n=A.size

	for i in range(n/2):
		print('frequency indices'),  -n/2 + 1+ i, n/2-i
		print('values'), A[i], A[-i-1]
		print('sum of imaginary part'), np.imag(A[i] + A[-i-1])
		x=A[i] + A[-i-1]
		print('should be zero'), np.imag(x)
	print('the zero frequency at '), n/2, '=', A[n/2],
	return 1

def J_k(k,P,param_matrix,nu=-2,P2=None,P_window=None, C_window=None,n_pad=500,verbose=False):

	# size of input array must be an even number
	if (k.size % 2 != 0):
		raise ValueError('Input array must contain an even number of elements.')

	alpha=param_matrix[:,0]
	beta=param_matrix[:,1]
	l_Bessel=param_matrix[:,2]
	type=param_matrix[:,3]

	N=k.size

	delta_L=(log(np.max(k))-log(np.min(k)))/(N-1)

	P_b=P*k**(-nu)

	if P_window is not None:
		# window the input power spectrum, so that at high and low k
		# the signal smoothly tappers to zero. This make the input
		# more "like" a periodic signal

		if (verbose):
			print('smoothing biased power spectrum')
		W=p_window(k,P_window[0],P_window[1])
		P_b=P_b*W

	if (n_pad !=None):
		# pad the edges. This helps with edge effects in Fourier space

		if (verbose):
			print('padding the input signal with'), n_pad, 'zeros.'
		id_pad=np.arange(k.size)
		k,P_b=pad_left(k,P_b,n_pad)
		_,P=pad_left(k,P,n_pad)
		k,P_b=pad_right(k,P_b,n_pad)
		_,P=pad_right(k,P,n_pad)
		N=k.size
		id_pad=id_pad+n_pad


	# I take the real Fourier transform and then take the conjugate to
	# obtain the negative frequencies.  The convolution latter in the code requires
	# the negative frequencies.
	c_m_positive=rfft(P_b)
	c_m_negative=np.conjugate(c_m_positive[1:])
	c_m=np.hstack((c_m_negative[::-1], c_m_positive))/float(N)

	# frequency integers
	n_c=c_m_positive.size
	m=np.arange(-n_c+1,n_c)


	# define eta_m and eta_n=eta_m
	omega=2*pi/(float(N)*delta_L)
	eta_m=omega*m

	# define l and tau_l
	n_l=c_m.size + c_m.size - 1
	l=l=np.arange(-n_l//2+1,n_l//2+1)
	tau_l=omega*l


	if (C_window != None):
		# window the Fourier coefficients.
		# This will damping the highest frequencies
		if (verbose):
			print('smoothing the Fourier coefficients')
		c_m=c_m*c_window(m,int(C_window*N//2.))

	# matrix for output
	A_out=np.zeros((param_matrix.shape[0],k.size))
	for i in range(param_matrix.shape[0]):

		sigma=l_Bessel[i]+1/2.

		# Define Q_m and Q_n and p
		# use eta_m for Q_n, the value is the same
		Q_m=3/2.+ nu + alpha[i] + 1j*eta_m
		Q_n=3/2.+ nu + beta[i] + 1j*eta_m
		p=-5-2*nu-alpha[i]-beta[i]

		g_m=g_m_vals(sigma,Q_m)

		if (type[i]==1):

			# this is the special case, Corresponding to the regularized version of
			# J_{2,-2,0,reg}(k)
			# get values for g_n
			# again use eta_m.
			s=2+nu + beta[i]
			Q_n=s+ 1j*eta_m

			g_n=gamsn(Q_n)

			#two_part_m=2**Q_m
			g_m=g_m*2.**Q_m

			# prefactor
			pf=(-1)**l_Bessel[i]/pi**3*np.sqrt(pi/2.)

			two_part_l=1

		else:
			g_n=g_m_vals(sigma,Q_n)
			# pre factor
			pf=(-1)**l_Bessel[i]/pi**2*2.**(2+2*nu+alpha[i]+beta[i])

			two_part_l=exp(1j*tau_l*log2)


		# convolve f_c and g_c
		C_l=np.convolve(c_m*g_m,c_m*g_n)

		# calculate h_l
		#arg=(p+1-1j*tau_l)
		h_l=gamsn(p+1-1j*tau_l)

		# multiply all l terms together
		C_l=C_l*h_l*two_part_l

		# set up to feed ifft an array ordered with l=0,1,...,-1,...,N/2-1
		c_plus=C_l[l>=0]
		c_minus=C_l[l< 0]

		C_l=np.hstack((c_plus[:-1],c_minus))
		A_k=ifft(C_l)*C_l.size # multiply by size to get rid of the normalization in ifft


		A_out[i,:]=np.real(A_k[::2])*pf*k**(-p-2) # note that you have to take every other element
												  # in A_k, due to the extended array created from the
												  # discrete convolution
	P_out=irfft(c_m[m>=0])*k**nu*float(N)
	if (n_pad !=0):
		# get rid of the elements created from padding
		P_out=P_out[id_pad]
		A_out=A_out[:,id_pad]

	return P_out, A_out
