import numpy as np
cimport numpy as cnp
import math
from mpmath import *

import scipy.special as sc
from scipy import linalg as sciLA
from scipy.linalg import eig


from libc.stdlib cimport rand, RAND_MAX
from libc.math cimport log, sqrt, cos, sin, acos, asin, exp
from libc.math cimport round

ctypedef cnp.float64_t DTYPE_t
DTYPE = np.float128
DTYPE2 = np.int32

mp.dps = 16; mp.pretty = False;

def rounding(double x):
	cdef double r
	r = round(x)
	return r
	
cdef double random_uniform():
	cdef double r = rand()
	return r / (RAND_MAX*1.0)

cdef double random_gaussian():
	cdef double u1, u2, w1
	u1 = random_uniform()
	u2 = random_uniform()
	if (u1==0.0): u1 = random_uniform()
	w1 = (-2.*log(u1))**0.5 * cos(2*3.14159265358979323846*u2)
#~ 	w2 = (-2.*log(u1))**0.5 * sin(2*3.14159265358979323846*u2)
	return w1

def myfunc ():
	print("Hello World")
	
def even_moments (double Pe, double alpha, double x0, double y0, double theta0, double dt, double time_duration):	
	# see main.py for a definition of various parameters
	cdef int t
	cdef double x,y,theta
	cdef int Nt
	cdef double factor
	cdef double r,r2
	
	cdef double even2,even4
	
	Nt = (int) (time_duration/dt)
	factor = sqrt(2*dt)
	
	x = x0
	y = y0
	theta = theta0
	even2 = 0.0
	even4 = 0.0
	for t in range(1,Nt+1):
		x += Pe * cos(theta) * dt - alpha*x*dt + factor *random_gaussian()
		y += Pe * sin(theta) * dt - alpha*y*dt + factor *random_gaussian()
		theta += factor * random_gaussian()
		
		r2 = x**2+y**2
		even2 += r2
		even4 += r2**2
	
	even2 = even2/(Nt)
	even4 = even4/(Nt)
		
	return even2,even4
			
def rcoschi_avg (double Pe, double alpha, double x0, double y0, double theta0, double dt, double time_duration):	
	# see main.py for a definition of various parameters
	cdef int t
	cdef double x,y,theta
	cdef int Nt
	cdef double factor
	cdef double r,phi
	
	cdef double even2,even4
	
	Nt = (int) (time_duration/dt)
	factor = sqrt(2*dt)
	
	x = x0
	y = y0
	theta = theta0
	rcoschi = 0.0
	for t in range(1,Nt+1):
		x += Pe * cos(theta) * dt - alpha*x*dt + factor *random_gaussian()
		y += Pe * sin(theta) * dt - alpha*y*dt + factor *random_gaussian()
		theta += factor * random_gaussian()
		
		r = sqrt(x**2+y**2)
		phi = acos(x/r)
		rcoschi += r * cos(theta-phi)
	
	rcoschi = rcoschi/(Nt)
		
	return rcoschi		
			
def spatial_distribution(double D_rot, double D, double v, double dt, double time_duration, double x0, double y0, double theta0, int Nparticles, double t1, double t2, double t3, double binning, int Nbins):		
	# see main.py for a definition of various parameters
	cdef int t,i,jx,jy
	cdef double x,y,theta
	cdef int Nt,Nt1,Nt2,Nt3
	cdef double factor,factor2,factor3,factor4
	cdef double r,r2
	
	Nt = (int) (time_duration/dt)
	Nt1 = (int) (t1/dt)
	Nt2 = (int) (t2/dt)
	Nt3 = (int) (t3/dt)
	
	factor = sqrt(2*D*dt)
	factor2 = sqrt(2*D_rot*dt)
	factor4 = v * dt
	
	
	data1 = np.zeros((Nbins,Nbins), dtype=DTYPE)
	data2 = np.zeros((Nbins,Nbins), dtype=DTYPE)
	data3 = np.zeros((Nbins,Nbins), dtype=DTYPE)
	
	
	for i in range(Nparticles):
		if (i % 100000 ==0): print i
		x = x0
		y = y0
		theta = theta0
		for t in range(1,Nt+1):
			x = x + factor4 * cos(theta) + factor *random_gaussian()
			y = y + factor4 * sin(theta) + factor *random_gaussian()
			theta = theta + factor2 *random_gaussian()
			
			
			if (x**2+y**2<=1):
				if (t == Nt1):
					jx = (int) (rounding(x/binning)) + Nbins/2
					jy = (int) (rounding(y/binning)) + Nbins/2
					if (jx >= 0 and jx<Nbins and jy >= 0 and jy<Nbins):
						data1[jx,jy] += 1.0
				elif (t == Nt2):
					jx = (int) (rounding(x/binning)) + Nbins/2
					jy = (int) (rounding(y/binning)) + Nbins/2
					if (jx >= 0 and jx<Nbins and jy >= 0 and jy<Nbins):
						data2[jx,jy] += 1.0
				elif (t == Nt3):
					jx = (int) (rounding(x/binning)) + Nbins/2
					jy = (int) (rounding(y/binning)) + Nbins/2
					if (jx >= 0 and jx<Nbins and jy >= 0 and jy<Nbins):
						data3[jx,jy] += 1.0
			else:
				break
	
	data1 = data1/Nparticles /binning**2
	data2 = data2/Nparticles /binning**2
	data3 = data3/Nparticles /binning**2
	
	return data1,data2,data3

def radial_distribution_avg(double D_rot, double D, double v, double dt, double time_duration, double r0, double phi0, int Nparticles, double t1, double t2, double t3, double binning, int Nbins):		
	# see main.py for a definition of various parameters
	cdef int t,i,jr
	cdef double x,y,theta
	cdef int Nt,Nt1,Nt2,Nt3
	cdef double factor,factor2,factor3,factor4
	cdef double r
	cdef double x0,y0
	
	Nt = (int) (time_duration/dt)
	Nt1 = (int) (t1/dt)
	Nt2 = (int) (t2/dt)
	Nt3 = (int) (t3/dt)
	
	factor = sqrt(2*D*dt)
	factor2 = sqrt(2*D_rot*dt)
	factor4 = v * dt
	
	rs = np.asarray([binning*i   for i in range(Nbins)])
	data1 = np.zeros((Nbins), dtype=DTYPE)
	data2 = np.zeros((Nbins), dtype=DTYPE)
	data3 = np.zeros((Nbins), dtype=DTYPE)
	
	x0 = r0*cos(phi0)
	y0 = r0*sin(phi0)
	
	
	for i in range(Nparticles):
		if (i % 100000 ==0): print i
		x = x0
		y = y0
		theta = random_uniform()*2*np.pi
		for t in range(1,Nt+1):
			x = x + factor4 * cos(theta) + factor *random_gaussian()
			y = y + factor4 * sin(theta) + factor *random_gaussian()
			theta = theta + factor2 *random_gaussian()
			
			r = sqrt(x**2+y**2)

			if (r<=1):
				if (t == Nt1):
					jr = (int) (rounding(r/binning))
					if (jr >= 0 and jr<Nbins):
						data1[jr] += 1.0
				elif (t == Nt2):
					jr = (int) (rounding(r/binning))
					if (jr >= 0 and jr<Nbins):
						data2[jr] += 1.0
				elif (t == Nt3):
					jr = (int) (rounding(r/binning))
					if (jr >= 0 and jr<Nbins):
						data3[jr] += 1.0
			else:
				break
	
	
	# leave it like that for normalization: survival probability = int_0^1 dr P(r)
	data1 = data1/Nparticles /binning
	data2 = data2/Nparticles /binning
	data3 = data3/Nparticles /binning
	#check
	suma1 = 0.0
	suma2 = 0.0
	suma3 = 0.0
	for i in range(len(rs)):
		suma1+=data1[i]*binning
		suma2+=data2[i]*binning
		suma3+=data3[i]*binning
	print 'Survival t1 = ',suma1
	print 'Survival t2 = ',suma2
	print 'Survival t3 = ',suma3
	
	# carefull in checking in the normalization: survival probability = int_0^1 dr 2 pi r P(r)
#~ 	for i in range(Nbins):
#~ 		if (i==0):
#~ 			r = rs[i+1]/4
#~ 			data1[i] = data1[i]/Nparticles /(binning*np.pi*r)
#~ 			data2[i] = data2[i]/Nparticles /(binning*np.pi*r)
#~ 			data3[i] = data3[i]/Nparticles /(binning*np.pi*r)
#~ 		else:
#~ 			r = rs[i]
#~ 			data1[i] = data1[i]/Nparticles /(binning*2*np.pi*r)
#~ 			data2[i] = data2[i]/Nparticles /(binning*2*np.pi*r)
#~ 			data3[i] = data3[i]/Nparticles /(binning*2*np.pi*r)
	#check
#~ 	suma1 = 0.0
#~ 	suma2 = 0.0
#~ 	suma3 = 0.0
#~ 	for i in range(len(rs)):
#~ 		suma1+=data1[i]*2*np.pi*rs[i]*binning
#~ 		suma2+=data2[i]*2*np.pi*rs[i]*binning
#~ 		suma3+=data3[i]*2*np.pi*rs[i]*binning
#~ 	print 'Survival t1 = ',suma1
#~ 	print 'Survival t2 = ',suma2
#~ 	print 'Survival t3 = ',suma3
	
	return data1,data2,data3

def survival(double D_rot, double D, double v, double dt, double r0, double phi0, double theta0, int Nparticles, double binning, int Nbins):		
	# see main.py for a definition of various parameters
	cdef int t,i,jt
	cdef double x,y,theta
	cdef int Nt
	cdef double factor,factor2,factor3,factor4
	cdef double r
	cdef double x0,y0
	cdef int Deltat
	
	ts = np.asarray([binning*i   for i in range(Nbins)])
	factor = sqrt(2*D*dt)
	factor2 = sqrt(2*D_rot*dt)
	factor4 = v * dt
	
	data1 = np.zeros((Nbins), dtype=DTYPE)
	
	x0 = r0*cos(phi0)
	y0 = r0*sin(phi0)
	
	Nt = (int) (ts[-1]/dt)
	Deltat = int (binning/dt)
	
	
	for i in range(Nparticles):
		if (i % 1000 ==0): print i
		x = x0
		y = y0
		theta = theta0
		data1[0]+=1.0
		
		for t in range(1,Nt+1):
			x = x + factor4 * cos(theta) + factor *random_gaussian()
			y = y + factor4 * sin(theta) + factor *random_gaussian()
			theta = theta + factor2 *random_gaussian()
			
			r = sqrt(x**2+y**2)

			if (r<=1):
				if (t%Deltat ==0):
					jt = (int) (t/Deltat)
					data1[jt]+=1.0
			else:
				break
	
	
	# leave it like that for normalization: survival probability = int_0^1 dr P(r)
	data1 = data1/Nparticles
	
	return data1

def fpt(double D_rot, double D, double v, double dt, double r0, double phi0, double theta0, int Nparticles, double binning, int Nbins):		
	# see main.py for a definition of various parameters
	cdef int t,i,jt
	cdef double x,y,theta
	cdef int Nt
	cdef double factor,factor2,factor3,factor4
	cdef double r
	cdef double x0,y0
	cdef int Deltat
	
	ts = np.asarray([binning*i   for i in range(Nbins)])
	factor = sqrt(2*D*dt)
	factor2 = sqrt(2*D_rot*dt)
	factor4 = v * dt
	
	data1 = np.zeros((Nbins), dtype=DTYPE)
	data2 = np.zeros((Nbins), dtype=DTYPE)
	
	x0 = r0*cos(phi0)
	y0 = r0*sin(phi0)
	
	Nt = (int) (ts[-1]/dt)
	Deltat = int (binning/dt)
	
	
	for i in range(Nparticles):
		if (i % 1000 ==0): print i
		x = x0
		y = y0
		theta = theta0
		data1[0]+=1.0
		
		for t in range(1,Nt+2):
			x = x + factor4 * cos(theta) + factor *random_gaussian()
			y = y + factor4 * sin(theta) + factor *random_gaussian()
			theta = theta + factor2 *random_gaussian()
			
			r = sqrt(x**2+y**2)

			if (r<=1):
				if (t%Deltat ==0):
					jt = (int) (t/Deltat)
					data1[jt]+=1.0
				if (t%Deltat ==1):
					jt = (int) (t/Deltat)
					data2[jt]+=1.0
			else:
				break
	
	
	# leave it like that for normalization: survival probability = int_0^1 dr P(r)
	data1 = data1/Nparticles
	data2 = data2/Nparticles
	
	data1 = (data1-data2)/dt
	
	return data1



########################################################################

def Jb(l,x):
	if (l>=0): return sc.jv(l,x)
	else: return (-1)**abs(l) * sc.jv(-l,x)

def jb(l,n):
	return sc.jn_zeros(l,n)[n-1]

def ev(n,l,j,Gamm):
	return -(jb(l,n)**2+Gamm*(j-l)**2)

def cp(n1,n,l):
	return jb(l,n)/(Jb(l+1,jb(l,n))) *jb(l+1,n1)*Jb(l+1,jb(l,n))*Jb(l,jb(l+1,n1))/(Jb(l+2,jb(l+1,n1))*(jb(l,n)**2-jb(l+1,n1)**2))

def cm(n1,n,l):
	return jb(l,n)/(Jb(l+1,jb(l,n))) *jb(l-1,n1)*Jb(l-1,jb(l,n))*Jb(l,jb(l-1,n1))/(Jb(l,jb(l-1,n1))*(jb(l,n)**2-jb(l-1,n1)**2))

def nn(a,Nn):
	return a%Nn+1

def ll(a,Nn,Nl):
	return int(a/Nn)-Nl

def aa(n,l,Nn,Nl):
	return (l+Nl)*Nn+(n-1)

def mat(a,b,j,Pe,Gamm,Nn,Nl):
	if a==b:
		return ev(nn(a,Nn),ll(a,Nn,Nl),j,Gamm)
	if ll(b,Nn,Nl)==(ll(a,Nn,Nl)+1):
		return Pe*cp(nn(b,Nn),nn(a,Nn),ll(a,Nn,Nl))
	if ll(b,Nn,Nl)==(ll(a,Nn,Nl)-1):
		return Pe*cm(nn(b,Nn),nn(a,Nn),ll(a,Nn,Nl))
	else:
		return 0.00000000
		
def psi(n,l,j,r,phi,theta):
#~ 	return np.exp(np.complex(0.0,l*phi)) * np.exp(np.complex(0.0,(j-l)*theta)) * Jb(l,jb(l,n)*r) / (sqrt(2.0)*np.pi * Jb(l+1,jb(l,n)))
	return Jb(l,jb(l,n)*r) / (sqrt(2.0)*np.pi * Jb(l+1,jb(l,n)))   * np.complex(cos(l*phi+(j-l)*theta),sin(l*phi+(j-l)*theta))
	
def psi_integrated(n,l,r,phi):
#~ 	return np.sqrt(2.0) * np.exp(np.complex(0.0,l*phi)) * Jb(l,jb(l,n)*r) / Jb(l+1,jb(l,n))
	return sqrt(2.0) * Jb(l,jb(l,n)*r) / Jb(l+1,jb(l,n)) * np.complex(cos(l*phi),sin(l*phi))
	
def psi_integrated_integrated(n,r):	
	return 2*sqrt(2.0) *np.pi * Jb(0,jb(0,n)*r) / Jb(1,jb(0,n))

def compare_eigenvalues(e1,e2):
	eps=0.0000001
	value = 0
	if (np.real(e1)<np.real(e2)+eps and np.real(e1)>np.real(e2)-eps and np.imag(e1)<np.imag(e2)+eps and np.imag(e1)>np.imag(e2)-eps):
		value=1
	return value

def spatial_distribution_numerics(double Gamm, double Pe, double r0, double phi0, double theta0, double t1, double t2,double t3, double binning, int Nbins):
	# size of matrices
	cdef int Nn=8
	cdef int Nl=7
	cdef int Ntot=Nn*(2*Nl+1)
	
	cdef int a,b,bp,i,j,k,n,l
	cdef double r,phi
	
	print('Size of matrix: ',Ntot)
	
	xs=np.asarray([ -(Nbins-1)/2*binning+ binning*i  for i in range(Nbins)])
	ys=np.asarray([ -(Nbins-1)/2*binning+ binning*i  for i in range(Nbins)])
	data1 = np.zeros((Nbins,Nbins),dtype=np.complex128)
	data2 = np.zeros((Nbins,Nbins),dtype=np.complex128)
	data3 = np.zeros((Nbins,Nbins),dtype=np.complex128)
	
	
	# different channels
	for j in range(-Nl,Nl+1):
		print j
		
		# initialize matrix, eigenvalues, eigenfunctions
		L = [[0.0 for i in range(Ntot)] for k in range(Ntot)]
		LT = [[0.0 for i in range(Ntot)] for k in range(Ntot)]
		for a in range(Ntot):
			for b in range(Ntot):
				L[b][a]=mat(a,b,j,Pe,Gamm,Nn,Nl)
				LT[a][b]=mat(a,b,j,Pe,Gamm,Nn,Nl)
		
		# find eigenvalues, and eigenvectors
		eigval, eigvec_right = sciLA.eig(L,left=False,right=True)
		eigval2, eigvec_left = sciLA.eig(LT,left=False,right=True)
		# correcly reorder left eigenvectors
		for a in range(Ntot):
			if (compare_eigenvalues(eigval[a],eigval2[a]) == 0):
				for b in range(Ntot):
					if (compare_eigenvalues(eigval[a],eigval2[b]) == 1):
						for c in range(Ntot):
							temp = eigvec_left[c,b]
							eigvec_left[c,b] = eigvec_left[c,a]
							eigvec_left[c,a] = temp
						temp = eigval2[b]
						eigval2[b] = eigval2[a]
						eigval2[a] = temp
						break
					if (b==Ntot-1): print 'problem with eigenvector decomposition'
		# correctly normalize eigenvectors
		M = np.dot(eigvec_left.T,eigvec_right)
		for a in range(Ntot):
#~ 			print M[a][a], eigval[a], eigval2[a],nn(a,Nn),ll(a,Nn,Nl)
			eigvec_left[:,a]=eigvec_left[:,a]/M[a][a]
#~ 		print ' '
		
		#################################
		factor1 = np.zeros((Ntot),dtype=np.complex128)
		for a in range(Ntot):
			for b in range(Ntot):
				n = nn(b,Nn)
				l = ll(b,Nn,Nl)
				factor1[a] += eigvec_left[b][a]*np.conj(psi(n,l,j,r0,phi0,theta0) )
		
		l = j
		for n in range(1,Nn+1):
			bp = aa(n,l,Nn,Nl)
			
			coeff1 = np.complex(0.0,0.0)
			coeff2 = np.complex(0.0,0.0)
			coeff3 = np.complex(0.0,0.0)
			
			for a in range(Ntot):
				coeff1 += np.exp(eigval[a] * t1) * eigvec_right[bp][a] * factor1[a]
				coeff2 += np.exp(eigval[a] * t2) * eigvec_right[bp][a] * factor1[a]
				coeff3 += np.exp(eigval[a] * t3) * eigvec_right[bp][a] * factor1[a]

			for i in range(len(xs)):
				x = xs[i]
				for k in range(len(ys)):
					y = ys[k]
			
					r = sqrt(x**2+y**2)
					phi = acos(x/r)
					if (y<0): phi = 2*np.pi - phi
					if (r==0.0): phi=0.0
					
					data1[i,k] += coeff1 * psi_integrated(n,l,r,phi)
					data2[i,k] += coeff2 * psi_integrated(n,l,r,phi)
					data3[i,k] += coeff3 * psi_integrated(n,l,r,phi)
		

	return np.real(data1), np.real(data2), np.real(data3)
	
def radial_distribution_avg_numerics(double Gamm, double Pe, double r0, double phi0, double t1, double t2,double t3, double binning, int Nbins):
	# size of matrices
	cdef int Nn=8
	cdef int Nl=7
	cdef int Ntot=Nn*(2*Nl+1)
	
	cdef int a,b1,b2,i,j,n2,n1
	cdef double r
	
	print('Size of matrix: ',Ntot)
	
	rs = np.asarray([binning*i   for i in range(Nbins)])
	data1 = np.zeros((Nbins),dtype=np.complex128)
	data2 = np.zeros((Nbins),dtype=np.complex128)
	data3 = np.zeros((Nbins),dtype=np.complex128)
	
	# only channel j=0
	j=0
		
	# initialize matrix, eigenvalues, eigenfunctions
	L = [[0.0 for i in range(Ntot)] for k in range(Ntot)]
	LT = [[0.0 for i in range(Ntot)] for k in range(Ntot)]
	for a in range(Ntot):
		for b in range(Ntot):
			L[b][a]=mat(a,b,j,Pe,Gamm,Nn,Nl)
			LT[a][b]=mat(a,b,j,Pe,Gamm,Nn,Nl)
	
	# find eigenvalues, and eigenvectors
	eigval, eigvec_right = sciLA.eig(L,left=False,right=True)
	eigval2, eigvec_left = sciLA.eig(LT,left=False,right=True)
	# correcly reorder left eigenvectors
	for a in range(Ntot):
		if (compare_eigenvalues(eigval[a],eigval2[a]) == 0):
			for b in range(Ntot):
				if (compare_eigenvalues(eigval[a],eigval2[b]) == 1):
					for c in range(Ntot):
						temp = eigvec_left[c,b]
						eigvec_left[c,b] = eigvec_left[c,a]
						eigvec_left[c,a] = temp
					temp = eigval2[b]
					eigval2[b] = eigval2[a]
					eigval2[a] = temp
					break
				if (b==Ntot-1): print 'problem with eigenvector decomposition'
	# correctly normalize eigenvectors
	M = np.dot(eigvec_left.T,eigvec_right)
	for a in range(Ntot):
#~ 			print M[a][a], eigval[a], eigval2[a],nn(a,Nn),ll(a,Nn,Nl)
		eigvec_left[:,a]=eigvec_left[:,a]/M[a][a]
#~ 		print ' '
	
	
	#################################
	factor1 = np.zeros((Ntot),dtype=np.complex128)
	for a in range(Ntot):
		for n1 in range(1,Nn+1):
			b1 = aa(n1,0,Nn,Nl)
			factor1[a] += eigvec_left[b1][a] * psi_integrated_integrated(n1,r0) /(4*np.pi*np.pi)
	
	
	for n2 in range(1,Nn+1):
		b2 = aa(n2,0,Nn,Nl)
		
		coeff1 = np.complex(0.0,0.0)
		coeff2 = np.complex(0.0,0.0)
		coeff3 = np.complex(0.0,0.0)
		
		for a in range(Ntot):
			coeff1 += np.exp(eigval[a] * t1) * eigvec_right[b2][a] * factor1[a]
			coeff2 += np.exp(eigval[a] * t2) * eigvec_right[b2][a] * factor1[a]
			coeff3 += np.exp(eigval[a] * t3) * eigvec_right[b2][a] * factor1[a]

		for i in range(len(rs)):
			r = rs[i]
			data1[i] += coeff1 * psi_integrated_integrated(n2,r)
			data2[i] += coeff2 * psi_integrated_integrated(n2,r)
			data3[i] += coeff3 * psi_integrated_integrated(n2,r)
	
	
	# carefull in checking in the normalization: survival probability = int_0^1 dr P(r)
	for i in range(len(rs)):
		data1[i]=data1[i]*rs[i]
		data2[i]=data2[i]*rs[i]
		data3[i]=data3[i]*rs[i]
	#check
	suma1 = 0.0
	suma2 = 0.0
	suma3 = 0.0
	for i in range(len(rs)):
		suma1+=data1[i]*binning
		suma2+=data2[i]*binning
		suma3+=data3[i]*binning
	print 'Survival t1 = ',suma1
	print 'Survival t2 = ',suma2
	print 'Survival t3 = ',suma3
	
	# carefull in checking in the normalization: survival probability = int_0^1 dr 2 pi r P(r)
#~ 	data1=data1/(2*np.pi)
#~ 	data2=data2/(2*np.pi)
#~ 	data3=data3/(2*np.pi)
#~ 	#check
#~ 	suma1 = 0.0
#~ 	suma2 = 0.0
#~ 	suma3 = 0.0
#~ 	for i in range(len(rs)):
#~ 		suma1+=data1[i]*2*np.pi*rs[i]*binning
#~ 		suma2+=data2[i]*2*np.pi*rs[i]*binning
#~ 		suma3+=data3[i]*2*np.pi*rs[i]*binning
#~ 	print 'Survival t1 = ',suma1
#~ 	print 'Survival t2 = ',suma2
#~ 	print 'Survival t3 = ',suma3	

	return np.real(data1), np.real(data2), np.real(data3)

def survival_numerics(double Gamm, double Pe, double r0, double phi0, double theta0, double binning, int Nbins):
	# size of matrices
	cdef int Nn=8
	cdef int Nl=7
	cdef int Ntot=Nn*(2*Nl+1)
	
	cdef int a,b1,b2,i,j,n2,n1,l1
	cdef double t
	
	print 'Numerics started!'
	print('Size of matrix: ',Ntot)
	
	ts = np.asarray([binning*i   for i in range(Nbins)])
	data1 = np.zeros((Nbins),dtype=np.complex128)
	
	# only channel j=0
	j=0
		
	# initialize matrix, eigenvalues, eigenfunctions
	L = [[0.0 for i in range(Ntot)] for k in range(Ntot)]
	LT = [[0.0 for i in range(Ntot)] for k in range(Ntot)]
	for a in range(Ntot):
		for b in range(Ntot):
			L[b][a]=mat(a,b,j,Pe,Gamm,Nn,Nl)
			LT[a][b]=mat(a,b,j,Pe,Gamm,Nn,Nl)
	
	# find eigenvalues, and eigenvectors
	eigval, eigvec_right = sciLA.eig(L,left=False,right=True)
	eigval2, eigvec_left = sciLA.eig(LT,left=False,right=True)
	# correcly reorder left eigenvectors
	for a in range(Ntot):
		if (compare_eigenvalues(eigval[a],eigval2[a]) == 0):
			for b in range(Ntot):
				if (compare_eigenvalues(eigval[a],eigval2[b]) == 1):
					for c in range(Ntot):
						temp = eigvec_left[c,b]
						eigvec_left[c,b] = eigvec_left[c,a]
						eigvec_left[c,a] = temp
					temp = eigval2[b]
					eigval2[b] = eigval2[a]
					eigval2[a] = temp
					break
				if (b==Ntot-1): print 'problem with eigenvector decomposition'
	# correctly normalize eigenvectors
	M = np.dot(eigvec_left.T,eigvec_right)
	for a in range(Ntot):
#~ 			print M[a][a], eigval[a], eigval2[a],nn(a,Nn),ll(a,Nn,Nl)
		eigvec_left[:,a]=eigvec_left[:,a]/M[a][a]
#~ 		print ' '
	
	
	#################################
	factor1 = np.zeros((Ntot),dtype=np.complex128)
	factor2 = np.zeros((Ntot),dtype=np.complex128)
	for a in range(Ntot):
		for n1 in range(1,Nn+1):
			for l1 in range(-Nl,Nl+1):
				b1 = aa(n1,l1,Nn,Nl)
				factor1[a] += eigvec_left[b1][a] * psi(n1,l1,0,r0,phi0,theta0)
	
	
		for n2 in range(1,Nn+1):
			b2 = aa(n2,0,Nn,Nl)
			factor2[a] += eigvec_right[b2][a]/jb(0,n2)
		
	for i in range(len(ts)):
		t = ts[i]
		for a in range(Ntot):
			data1[i] += np.exp(eigval[a] * t) * factor1[a] * factor2[a]
		data1[i] *= 2 * np.pi * np.sqrt(2.0)
		
	if (ts[0]==0.0):
		data1[0] = np.complex(1.0,0.0)
	
	print 'Numerics finished!'
	
	return np.real(data1)

def fpt_numerics(double Gamm, double Pe, double r0, double phi0, double theta0, double binning, int Nbins):
	# size of matrices
	cdef int Nn=10
	cdef int Nl=9
	cdef int Ntot=Nn*(2*Nl+1)
	
	cdef int a,b1,b2,i,j,n2,n1,l1
	cdef double t
	
	print 'Numerics started!'
	print('Size of matrix: ',Ntot)
	
	ts = np.asarray([binning*i   for i in range(Nbins)])
	data1 = np.zeros((Nbins),dtype=np.complex128)
	
	# only channel j=0
	j=0
		
	# initialize matrix, eigenvalues, eigenfunctions
	L = [[0.0 for i in range(Ntot)] for k in range(Ntot)]
	LT = [[0.0 for i in range(Ntot)] for k in range(Ntot)]
	for a in range(Ntot):
		for b in range(Ntot):
			L[b][a]=mat(a,b,j,Pe,Gamm,Nn,Nl)
			LT[a][b]=mat(a,b,j,Pe,Gamm,Nn,Nl)
	
	# find eigenvalues, and eigenvectors
	eigval, eigvec_right = sciLA.eig(L,left=False,right=True)
	eigval2, eigvec_left = sciLA.eig(LT,left=False,right=True)
	# correcly reorder left eigenvectors
	for a in range(Ntot):
		if (compare_eigenvalues(eigval[a],eigval2[a]) == 0):
			for b in range(Ntot):
				if (compare_eigenvalues(eigval[a],eigval2[b]) == 1):
					for c in range(Ntot):
						temp = eigvec_left[c,b]
						eigvec_left[c,b] = eigvec_left[c,a]
						eigvec_left[c,a] = temp
					temp = eigval2[b]
					eigval2[b] = eigval2[a]
					eigval2[a] = temp
					break
				if (b==Ntot-1): print 'problem with eigenvector decomposition'
	# correctly normalize eigenvectors
	M = np.dot(eigvec_left.T,eigvec_right)
	for a in range(Ntot):
#~ 			print M[a][a], eigval[a], eigval2[a],nn(a,Nn),ll(a,Nn,Nl)
		eigvec_left[:,a]=eigvec_left[:,a]/M[a][a]
#~ 		print ' '
	
	
	#################################
	factor1 = np.zeros((Ntot),dtype=np.complex128)
	factor2 = np.zeros((Ntot),dtype=np.complex128)
	for a in range(Ntot):
		for n1 in range(1,Nn+1):
			for l1 in range(-Nl,Nl+1):
				b1 = aa(n1,l1,Nn,Nl)
				factor1[a] += eigvec_left[b1][a] * psi(n1,l1,0,r0,phi0,theta0)
	
	
		for n2 in range(1,Nn+1):
			b2 = aa(n2,0,Nn,Nl)
			factor2[a] += eigvec_right[b2][a]/jb(0,n2)
		
	for i in range(len(ts)):
		t = ts[i]
		for a in range(Ntot):
			data1[i] -= eigval[a] * np.exp(eigval[a] * t) * factor1[a] * factor2[a]
		data1[i] *= 2 * np.pi * np.sqrt(2.0)
		
	if (ts[0]==0.0):
		data1[0] = np.complex(0.0,0.0)
	
	print 'Numerics finished!'
	
	return np.real(data1)

def indexn(n):
	return n-1

def indexl(l,Nl):
	return l+Nl

def fq1(n,l,j,n1,l1,j1,t,Gamm):
	if (ev(n,l,j,Gamm) == ev(n1,l1,j1,Gamm) ):
		value = t*np.exp(ev(n,l,j,Gamm)*t)
	else:
		value = (np.exp(ev(n1,l1,j1,Gamm)*t) - np.exp(ev(n,l,j,Gamm)*t) ) / (ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))
	return value
	
def fq1_tder(n,l,j,n1,l1,j1,t,Gamm):
	if (ev(n,l,j,Gamm) == ev(n1,l1,j1,Gamm) ):
		value = np.exp(ev(n,l,j,Gamm)*t) + ev(n,l,j,Gamm)*t*np.exp(ev(n,l,j,Gamm)*t)
	else:
		value = (ev(n1,l1,j1,Gamm)*np.exp(ev(n1,l1,j1,Gamm)*t) - ev(n,l,j,Gamm)*np.exp(ev(n,l,j,Gamm)*t) ) / (ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))
	return value	
	
def fq2(n,l,j,n1,l1,j1,n2,l2,j2,t,Gamm):
	if (ev(n1,l1,j1,Gamm) == ev(n2,l2,j2,Gamm) ):
		if (ev(n,l,j,Gamm) == ev(n1,l1,j1,Gamm)):
			value = t**2/2
		else:
			value = (1+np.exp((ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))*t) * ((ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))*t-1))/(ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))**2
	else:
		if (  ev(n,l,j,Gamm) == ev(n1,l1,j1,Gamm) ):
			value = (np.exp((ev(n2,l2,j2,Gamm)-ev(n,l,j,Gamm))*t) -1 -(ev(n2,l2,j2,Gamm)-ev(n,l,j,Gamm))*t)/(ev(n2,l2,j2,Gamm)-ev(n,l,j,Gamm))**2
		elif ( ev(n,l,j,Gamm) == ev(n2,l2,j2,Gamm)  ):
			value = (np.exp((ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))*t) -1 -(ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))*t)/(ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))**2
		else:
			value = (np.exp((ev(n2,l2,j2,Gamm)-ev(n,l,j,Gamm))*t) -1)/(ev(n2,l2,j2,Gamm)-ev(n,l,j,Gamm)) - (np.exp((ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))*t) -1)/(ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))
			value = value/(ev(n2,l2,j2,Gamm)-ev(n1,l1,j1,Gamm))
	return value

def fq2_tder(n,l,j,n1,l1,j1,n2,l2,j2,t,Gamm):
	if (ev(n1,l1,j1,Gamm) == ev(n2,l2,j2,Gamm) ):
		if (ev(n,l,j,Gamm) == ev(n1,l1,j1,Gamm)):
			value = t
		else:
			value = (np.exp((ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))*t) * ((ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))*t-1) + np.exp((ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))*t))/(ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))
	else:
		if (  ev(n,l,j,Gamm) == ev(n1,l1,j1,Gamm) ):
			value = (np.exp((ev(n2,l2,j2,Gamm)-ev(n,l,j,Gamm))*t) - 1)/(ev(n2,l2,j2,Gamm)-ev(n,l,j,Gamm))
		elif ( ev(n,l,j,Gamm) == ev(n2,l2,j2,Gamm)  ):
			value = (np.exp((ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))*t) -1)/(ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))
		else:
			value = np.exp((ev(n2,l2,j2,Gamm)-ev(n,l,j,Gamm))*t) - np.exp((ev(n1,l1,j1,Gamm)-ev(n,l,j,Gamm))*t)
			value = value/(ev(n2,l2,j2,Gamm)-ev(n1,l1,j1,Gamm))
	return value

	
def spatial_distribution_analytics (double Gamm, double Pe, double r0, double phi0, double theta0, double t1, double t2,double t3, double binning, int Nbins):
	cdef int Nn=8
	cdef int Nl=7
	cdef int qmax=2
	
	cdef int i,j,k,n,l,q,n1,n2,l1,l2
	
	print 'analytics starting!'
	
	xs=np.asarray([ -(Nbins-1)/2*binning+ binning*i  for i in range(Nbins)])
	ys=np.asarray([ -(Nbins-1)/2*binning+ binning*i  for i in range(Nbins)])
	data1 = np.zeros((Nbins,Nbins),dtype=np.complex128)
	data2 = np.zeros((Nbins,Nbins),dtype=np.complex128)
	data3 = np.zeros((Nbins,Nbins),dtype=np.complex128)
	
	
	M1q = np.zeros((Nn,2*Nl+1,2*Nl+1,qmax+1),dtype=np.complex128)
	M2q = np.zeros((Nn,2*Nl+1,2*Nl+1,qmax+1),dtype=np.complex128)
	M3q = np.zeros((Nn,2*Nl+1,2*Nl+1,qmax+1),dtype=np.complex128)
	M1 = np.zeros((Nn,2*Nl+1,2*Nl+1),dtype=np.complex128)
	M2 = np.zeros((Nn,2*Nl+1,2*Nl+1),dtype=np.complex128)
	M3 = np.zeros((Nn,2*Nl+1,2*Nl+1),dtype=np.complex128)
	
	
	for q in range(qmax+1):
		if (q==0):
			for n in range(1,Nn+1):
				print 'q=0 n=',n
				for l in range(-Nl,Nl+1):
					j=l
					M1q[indexn(n),indexl(l,Nl),indexl(j,Nl),0] = np.exp(ev(n,l,j,Gamm)*t1) * np.conj(psi(n,l,j,r0,phi0,theta0))
					M2q[indexn(n),indexl(l,Nl),indexl(j,Nl),0] = np.exp(ev(n,l,j,Gamm)*t2) * np.conj(psi(n,l,j,r0,phi0,theta0))
					M3q[indexn(n),indexl(l,Nl),indexl(j,Nl),0] = np.exp(ev(n,l,j,Gamm)*t3) * np.conj(psi(n,l,j,r0,phi0,theta0))
						
		elif (q==1):
			for n in range(1,Nn+1):
				print 'q=1 n=',n
				for l in range(-Nl,Nl+1):
					j=l
						
					for n1 in range(1,Nn+1):
						factor = cp(n,n1,l-1) * np.conj(psi(n1,l-1,j,r0,phi0,theta0))
						M1q[indexn(n),indexl(l,Nl),indexl(j,Nl),1] += factor * fq1(n,l,j,n1,l-1,j,t1,Gamm)
						M2q[indexn(n),indexl(l,Nl),indexl(j,Nl),1] += factor * fq1(n,l,j,n1,l-1,j,t2,Gamm)
						M3q[indexn(n),indexl(l,Nl),indexl(j,Nl),1] += factor * fq1(n,l,j,n1,l-1,j,t3,Gamm)
						
						factor = cm(n,n1,l+1) * np.conj(psi(n1,l+1,j,r0,phi0,theta0))
						M1q[indexn(n),indexl(l,Nl),indexl(j,Nl),1] += factor * fq1(n,l,j,n1,l+1,j,t1,Gamm)
						M2q[indexn(n),indexl(l,Nl),indexl(j,Nl),1] += factor * fq1(n,l,j,n1,l+1,j,t2,Gamm)
						M3q[indexn(n),indexl(l,Nl),indexl(j,Nl),1] += factor * fq1(n,l,j,n1,l+1,j,t3,Gamm)
			
		elif (q==2):
			for n in range(1,Nn+1):
				print 'q=2 n=',n
				for l in range(-Nl,Nl+1):
					j=l
					for n1 in range(1,Nn+1):
						for n2 in range(1,Nn+1):
							factor = cp(n,n1,l-1) * cp(n1,n2,l-2) * np.conj(psi(n2,l-2,j,r0,phi0,theta0))
							M1q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l-2,j,t1,Gamm)
							M2q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l-2,j,t2,Gamm)
							M3q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l-2,j,t3,Gamm)
							
							factor = cp(n,n1,l-1) * cm(n1,n2,l) * np.conj(psi(n2,l,j,r0,phi0,theta0))
							M1q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l,j,t1,Gamm)
							M2q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l,j,t2,Gamm)
							M3q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l,j,t3,Gamm)
							
							factor = cm(n,n1,l+1) * cp(n1,n2,l) * np.conj(psi(n2,l,j,r0,phi0,theta0))
							M1q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l,j,t1,Gamm)
							M2q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l,j,t2,Gamm)
							M3q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l,j,t3,Gamm)
							
							factor = cm(n,n1,l+1) * cm(n1,n2,l+2) * np.conj(psi(n2,l+2,j,r0,phi0,theta0))
							M1q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l+2,j,t1,Gamm)
							M2q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l+2,j,t2,Gamm)
							M3q[indexn(n),indexl(l,Nl),indexl(j,Nl),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l+2,j,t3,Gamm)
					M1q[indexn(n),indexl(l,Nl),indexl(j,Nl),2] *= np.exp(ev(n,l,j,Gamm)*t1)
					M2q[indexn(n),indexl(l,Nl),indexl(j,Nl),2] *= np.exp(ev(n,l,j,Gamm)*t2)
					M3q[indexn(n),indexl(l,Nl),indexl(j,Nl),2] *= np.exp(ev(n,l,j,Gamm)*t3)
			
		else:
			print 'not ready yet for this qmax = ',qmax
	
		
		
		
	
	
	for n in range(1,Nn+1):
			for l in range(-Nl,Nl+1):
				for j in range(-Nl,Nl+1):
					for q in range(qmax+1):
						M1[indexn(n),indexl(l,Nl),indexl(j,Nl)] += Pe**q * M1q[indexn(n),indexl(l,Nl),indexl(j,Nl),q]
						M2[indexn(n),indexl(l,Nl),indexl(j,Nl)] += Pe**q * M2q[indexn(n),indexl(l,Nl),indexl(j,Nl),q]
						M3[indexn(n),indexl(l,Nl),indexl(j,Nl)] += Pe**q * M3q[indexn(n),indexl(l,Nl),indexl(j,Nl),q]
						
						
	
	print 'computation of Ms finished!'
	
	
	for i in range(len(xs)):
		x = xs[i]
		for k in range(len(ys)):
			y = ys[k]
			
			r = sqrt(x**2+y**2)
			
			if (r==0.0): 
				phi=0.0
			else:
				phi = acos(x/r)
				if (y<0): phi = 2*np.pi - phi
				
			if (r>=1):
				data1[i,k]=np.complex(0.0,0.0)
				data2[i,k]=np.complex(0.0,0.0)
				data3[i,k]=np.complex(0.0,0.0)
			else:
				if (y==0.0): print x
				for n in range(1,Nn+1):
					for l in range(-Nl,Nl+1):
						data1[i,k] += M1[indexn(n),indexl(l,Nl),indexl(l,Nl)] * psi_integrated(n,l,r,phi)
						data2[i,k] += M2[indexn(n),indexl(l,Nl),indexl(l,Nl)] * psi_integrated(n,l,r,phi)
						data3[i,k] += M3[indexn(n),indexl(l,Nl),indexl(l,Nl)] * psi_integrated(n,l,r,phi)
					
	print 'analytics finished!'
					
	
	
	return np.real(data1), np.real(data2), np.real(data3)
	
def radial_distribution_avg_analytics (double Gamm, double Pe, double r0, double phi0, double t1, double t2,double t3, double binning, int Nbins):
	# size of matrices
	cdef int Nn=8
	cdef int Nl=7
	cdef int qmax=2
	
	cdef double r
	cdef int i,j,k,n,l,q,n1,n2,l1,l2
	
	rs = np.asarray([binning*i   for i in range(Nbins)])
	data1 = np.zeros((Nbins),dtype=np.complex128)
	data2 = np.zeros((Nbins),dtype=np.complex128)
	data3 = np.zeros((Nbins),dtype=np.complex128)
	
	print 'analytics starting!'
	
	rs = np.asarray([binning*i   for i in range(Nbins)])
	data1 = np.zeros((Nbins),dtype=np.complex128)
	data2 = np.zeros((Nbins),dtype=np.complex128)
	data3 = np.zeros((Nbins),dtype=np.complex128)
	
	
	M1q = np.zeros((Nn,2*Nl+1,qmax+1),dtype=np.complex128)
	M2q = np.zeros((Nn,2*Nl+1,qmax+1),dtype=np.complex128)
	M3q = np.zeros((Nn,2*Nl+1,qmax+1),dtype=np.complex128)
	M1 = np.zeros((Nn,2*Nl+1),dtype=np.complex128)
	M2 = np.zeros((Nn,2*Nl+1),dtype=np.complex128)
	M3 = np.zeros((Nn,2*Nl+1),dtype=np.complex128)
	
	j=0
	
	for q in range(qmax+1):
		if (q==0):
			for n in range(1,Nn+1):
				print 'q=0 n=',n
				for l in range(-Nl,Nl+1):
					M1q[indexn(n),indexl(l,Nl),0] = np.exp(ev(n,0,0,Gamm)*t1) * psi_integrated_integrated(n,r0) /(4*np.pi**2)
					M2q[indexn(n),indexl(l,Nl),0] = np.exp(ev(n,0,0,Gamm)*t2) * psi_integrated_integrated(n,r0) /(4*np.pi**2)
					M3q[indexn(n),indexl(l,Nl),0] = np.exp(ev(n,0,0,Gamm)*t3) * psi_integrated_integrated(n,r0) /(4*np.pi**2)
						
		elif (q==1):
			print 'q=1 No contribution to radial probability from this order'
				
			
		elif (q==2):
			for n in range(1,Nn+1):
				print 'q=2 n=',n
				l=0
				for n1 in range(1,Nn+1):
					for n2 in range(1,Nn+1):
						factor = cp(n,n1,l-1) * cm(n1,n2,l) * psi_integrated_integrated(n2,r0) /(4*np.pi**2)
						M1q[indexn(n),indexl(l,Nl),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l,j,t1,Gamm)
						M2q[indexn(n),indexl(l,Nl),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l,j,t2,Gamm)
						M3q[indexn(n),indexl(l,Nl),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l,j,t3,Gamm)
						factor = cm(n,n1,l+1) * cp(n1,n2,l) * psi_integrated_integrated(n2,r0) /(4*np.pi**2)
						M1q[indexn(n),indexl(l,Nl),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l,j,t1,Gamm)
						M2q[indexn(n),indexl(l,Nl),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l,j,t2,Gamm)
						M3q[indexn(n),indexl(l,Nl),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l,j,t3,Gamm)
					

				M1q[indexn(n),indexl(l,Nl),2] *= np.exp(ev(n,l,j,Gamm)*t1)
				M2q[indexn(n),indexl(l,Nl),2] *= np.exp(ev(n,l,j,Gamm)*t2)
				M3q[indexn(n),indexl(l,Nl),2] *= np.exp(ev(n,l,j,Gamm)*t3)
		
		else:
			print 'not ready yet for this qmax = ',qmax
	

	for n in range(1,Nn+1):
			for l in range(-Nl,Nl+1):
				for q in range(qmax+1):
					M1[indexn(n),indexl(l,Nl)] += Pe**q * M1q[indexn(n),indexl(l,Nl),q]
					M2[indexn(n),indexl(l,Nl)] += Pe**q * M2q[indexn(n),indexl(l,Nl),q]
					M3[indexn(n),indexl(l,Nl)] += Pe**q * M3q[indexn(n),indexl(l,Nl),q]
	
	print 'computation of Ms finished!'
	
	
	for i in range(len(rs)):
		r = rs[i]
		for n in range(1,Nn+1):
			data1[i] += M1[indexn(n),indexl(0,Nl)] * psi_integrated_integrated(n,r)
			data2[i] += M2[indexn(n),indexl(0,Nl)] * psi_integrated_integrated(n,r)
			data3[i] += M3[indexn(n),indexl(0,Nl)] * psi_integrated_integrated(n,r)
	
	

	# carefull in checking in the normalization: survival probability = int_0^1 dr P(r)
	for i in range(len(rs)):
		data1[i]=data1[i]*rs[i]
		data2[i]=data2[i]*rs[i]
		data3[i]=data3[i]*rs[i]
	#check
	suma1 = 0.0
	suma2 = 0.0
	suma3 = 0.0
	for i in range(len(rs)):
		suma1+=data1[i]*binning
		suma2+=data2[i]*binning
		suma3+=data3[i]*binning
	print 'Survival t1 = ',suma1
	print 'Survival t2 = ',suma2
	print 'Survival t3 = ',suma3
	
	# carefull in checking in the normalization: survival probability = int_0^1 dr 2 pi r P(r)
#~ 	data1=data1/(2*np.pi)
#~ 	data2=data2/(2*np.pi)
#~ 	data3=data3/(2*np.pi)
#~ 	#check
#~ 	suma1 = 0.0
#~ 	suma2 = 0.0
#~ 	suma3 = 0.0
#~ 	for i in range(len(rs)):
#~ 		suma1+=data1[i]*2*np.pi*rs[i]*binning
#~ 		suma2+=data2[i]*2*np.pi*rs[i]*binning
#~ 		suma3+=data3[i]*2*np.pi*rs[i]*binning
#~ 	print 'Survival t1 = ',suma1
#~ 	print 'Survival t2 = ',suma2
#~ 	print 'Survival t3 = ',suma3
	
	
	
	
	
	
	
	
					
	print 'analytics finished!'
					
	
	
	return np.real(data1), np.real(data2), np.real(data3)
	
def survival_analytics (double Gamm, double Pe, double r0, double phi0, double theta0, double binning, int Nbins):
	# size of matrices
	cdef int Nn=8
	cdef int Nl=7
	cdef int qmax=2
	
	cdef int i,j,k,n,l,q,n1,n2,l1,l2
	
	ts = np.asarray([binning*i   for i in range(Nbins)])
	data1 = np.zeros((Nbins),dtype=np.complex128)
	
	print 'analytics starting!'
	
	for i in range(len(ts)):
		t1 = ts[i]
		print t1
	
		M1q = np.zeros((Nn,qmax+1),dtype=np.complex128)
		M1 = np.zeros((Nn),dtype=np.complex128)
		
		j=0
		l=0
		
		for q in range(qmax+1):
			if (q==0):
				for n in range(1,Nn+1):
					M1q[indexn(n),0] = np.exp(ev(n,l,j,Gamm)*t1) * np.conj(psi(n,l,j,r0,phi0,theta0))
							
			elif (q==1):
				for n in range(1,Nn+1):
					for n1 in range(1,Nn+1):
						factor = cp(n,n1,l-1) * np.conj(psi(n1,l-1,j,r0,phi0,theta0))
						M1q[indexn(n),1] += factor * fq1(n,l,j,n1,l-1,j,t1,Gamm)
						
						factor = cm(n,n1,l+1) * np.conj(psi(n1,l+1,j,r0,phi0,theta0))
						M1q[indexn(n),1] += factor * fq1(n,l,j,n1,l+1,j,t1,Gamm)
				
			elif (q==2):
				for n in range(1,Nn+1):
					for n1 in range(1,Nn+1):
						for n2 in range(1,Nn+1):
							factor = cp(n,n1,l-1) * cp(n1,n2,l-2) * np.conj(psi(n2,l-2,j,r0,phi0,theta0))
							M1q[indexn(n),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l-2,j,t1,Gamm)
							
							factor = cp(n,n1,l-1) * cm(n1,n2,l) * np.conj(psi(n2,l,j,r0,phi0,theta0))
							M1q[indexn(n),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l,j,t1,Gamm)
							
							factor = cm(n,n1,l+1) * cp(n1,n2,l) * np.conj(psi(n2,l,j,r0,phi0,theta0))
							M1q[indexn(n),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l,j,t1,Gamm)
							
							factor = cm(n,n1,l+1) * cm(n1,n2,l+2) * np.conj(psi(n2,l+2,j,r0,phi0,theta0))
							M1q[indexn(n),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l+2,j,t1,Gamm)
					M1q[indexn(n),2] *= np.exp(ev(n,l,j,Gamm)*t1)
				
			else:
				print 'not ready yet for this qmax = ',qmax
		

		for n in range(1,Nn+1):
			for q in range(qmax+1):
				M1[indexn(n)] += Pe**q * M1q[indexn(n),q]
		
		
		for n in range(1,Nn+1):
			data1[i] += 2*np.sqrt(2)*np.pi*M1[indexn(n)] / jb(0,n)
	
	
#~ 		if (i==1 or ts[i]==0.1 or ts[i]==0.2 or ts[i]==0.3 or ts[i]==0.4):
#~ 			print 'check at t = ',ts[i]
#~ 			print np.real(2*np.sqrt(2)*np.pi*M1[indexn(1)] / jb(0,1)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(2)] / jb(0,2)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(3)] / jb(0,3)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(4)] / jb(0,4)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(5)] / jb(0,5)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(6)] / jb(0,6)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(7)] / jb(0,7)/data1[i])
	
	if (ts[0]==0.0):
		data1[0] = np.complex(1.0,0.0)
		
	print 'analytics finished!'
					
	
	
	return np.real(data1)
	
def fpt_analytics (double Gamm, double Pe, double r0, double phi0, double theta0, double binning, int Nbins):
	# size of matrices
	cdef int Nn=8
	cdef int Nl=7
	cdef int qmax=2
	
	cdef int i,j,k,n,l,q,n1,n2,l1,l2
	
	ts = np.asarray([binning*i   for i in range(Nbins)])
	data1 = np.zeros((Nbins),dtype=np.complex128)
	
	print 'fpt analytics starting!'
	
	for i in range(len(ts)):
		t1 = ts[i]
		print t1
	
		M1q = np.zeros((Nn,qmax+1),dtype=np.complex128)
		M1 = np.zeros((Nn),dtype=np.complex128)
		
		partialM1q = np.zeros((Nn,qmax+1),dtype=np.complex128)
		
		j=0
		l=0
		
		for q in range(qmax+1):
			if (q==0):
				for n in range(1,Nn+1):
					M1q[indexn(n),0] = ev(n,l,j,Gamm) * np.exp(ev(n,l,j,Gamm)*t1) * np.conj(psi(n,l,j,r0,phi0,theta0))
							
			elif (q==1):
				for n in range(1,Nn+1):
					for n1 in range(1,Nn+1):
						factor = cp(n,n1,l-1) * np.conj(psi(n1,l-1,j,r0,phi0,theta0))
						M1q[indexn(n),1] += factor * fq1_tder(n,l,j,n1,l-1,j,t1,Gamm)
						
						factor = cm(n,n1,l+1) * np.conj(psi(n1,l+1,j,r0,phi0,theta0))
						M1q[indexn(n),1] += factor * fq1_tder(n,l,j,n1,l+1,j,t1,Gamm)
				
			elif (q==2):
				for n in range(1,Nn+1):
					for n1 in range(1,Nn+1):
						for n2 in range(1,Nn+1):
							factor = cp(n,n1,l-1) * cp(n1,n2,l-2) * np.conj(psi(n2,l-2,j,r0,phi0,theta0))
							M1q[indexn(n),2]+= factor * fq2_tder(n,l,j,n1,l-1,j,n2,l-2,j,t1,Gamm)
							partialM1q[indexn(n),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l-2,j,t1,Gamm)
							
							factor = cp(n,n1,l-1) * cm(n1,n2,l) * np.conj(psi(n2,l,j,r0,phi0,theta0))
							M1q[indexn(n),2]+= factor * fq2_tder(n,l,j,n1,l-1,j,n2,l,j,t1,Gamm)
							partialM1q[indexn(n),2]+= factor * fq2(n,l,j,n1,l-1,j,n2,l,j,t1,Gamm)
							
							factor = cm(n,n1,l+1) * cp(n1,n2,l) * np.conj(psi(n2,l,j,r0,phi0,theta0))
							M1q[indexn(n),2]+= factor * fq2_tder(n,l,j,n1,l+1,j,n2,l,j,t1,Gamm)
							partialM1q[indexn(n),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l,j,t1,Gamm)
							
							factor = cm(n,n1,l+1) * cm(n1,n2,l+2) * np.conj(psi(n2,l+2,j,r0,phi0,theta0))
							M1q[indexn(n),2]+= factor * fq2_tder(n,l,j,n1,l+1,j,n2,l+2,j,t1,Gamm)
							partialM1q[indexn(n),2]+= factor * fq2(n,l,j,n1,l+1,j,n2,l+2,j,t1,Gamm)
					M1q[indexn(n),2] = np.exp(ev(n,l,j,Gamm)*t1)*M1q[indexn(n),2] + ev(n,l,j,Gamm) * np.exp(ev(n,l,j,Gamm)*t1) * partialM1q[indexn(n),2]
			else:
				print 'not ready yet for this qmax = ',qmax
		

		for n in range(1,Nn+1):
			for q in range(qmax+1):
				M1[indexn(n)] += Pe**q * M1q[indexn(n),q]
		
		
		for n in range(1,Nn+1):
			data1[i] -= 2*np.sqrt(2)*np.pi*M1[indexn(n)] / jb(0,n)
	
	
#~ 		if (i==1 or ts[i]==0.1 or ts[i]==0.2 or ts[i]==0.3 or ts[i]==0.4):
#~ 			print 'check at t = ',ts[i]
#~ 			print np.real(2*np.sqrt(2)*np.pi*M1[indexn(1)] / jb(0,1)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(2)] / jb(0,2)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(3)] / jb(0,3)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(4)] / jb(0,4)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(5)] / jb(0,5)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(6)] / jb(0,6)/data1[i]),np.real(2*np.sqrt(2)*np.pi*M1[indexn(7)] / jb(0,7)/data1[i])
	
	if (ts[0]==0.0):
		data1[0] = np.complex(0.0,0.0)
		
	print 'fpt analytics finished!'
					
	
	
	return np.real(data1)