import numpy as np
cimport numpy as cnp
import math
from mpmath import *
import scipy.special
from scipy import linalg as sciLA

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 k, double mu, double dt, double time_duration, double x0, double y0, double theta0, int Nstories, 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)
	factor3 = mu * k * 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(Nstories):
		if (i % 100 ==0): print i
		x = x0
		y = y0
		theta = theta0
		for t in range(1,Nt+1):
			x = x + factor4 * cos(theta) - factor3 * x + factor *random_gaussian()
			y = y + factor4 * sin(theta) - factor3 * y + factor *random_gaussian()
			theta = theta + factor2 *random_gaussian()
			
			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
	
	data1 = data1/Nstories /binning**2
	data2 = data2/Nstories /binning**2
	data3 = data3/Nstories /binning**2
	
	return data1,data2,data3
	
def landa (double D_rot, double mu, double k, int n, int l, int m):
	return D_rot*(m-l)**2 + mu*k*(2*n+abs(l))
	
def function1 (int n, int l):
	cdef double value
	value =1.0
	for i in range(n+1,n+l+1):
		value = value/(i*1.)
	return sqrt(value)
	
def psi (double r, double phi, double theta, double d, int n, int l, int m):
	cdef double prefactor
	prefactor = d**(0.5*abs(l)) * function1(n,abs(l)) * r**(abs(l)) * scipy.special.eval_genlaguerre(n, abs(l), d*r**2)
	return  prefactor * cos(l*phi-l*theta+m*theta), prefactor * sin(l*phi-l*theta+m*theta)
	
def psi_integrated_theta (double r, double phi, double d, int n, int l):
	# integreting over theta select m=l
	cdef double prefactor
	prefactor = 2* np.pi* d**(0.5*abs(l)) * function1(n,abs(l)) * r**(abs(l)) * scipy.special.eval_genlaguerre(n, abs(l), d*r**2)
	return  prefactor * cos(l*phi), prefactor * sin(l*phi)
	
def peq (double r, double d):
	return  d*exp(-d*r**2)/(2*np.pi**2)

def spatial_distribution_numerics_v1(double D_rot, double D, double v, double k, double mu, double x0, double y0, double theta0, double t, double binning, int Nbins, int qmax):		
	# see main.py for a definition of various parameters
	cdef int n,l,m
	cdef double x,y,theta
	cdef double r,phi
	cdef double d = mu*k/(2*D)
	cdef double a = v*sqrt(d)
	cdef int n_max = qmax/2
	cdef int l_max=2*n_max
	cdef int m_max = l_max
#~ 	cdef int m_max = 0
	# m_max should be equal to l_max
	cdef int q,p,s
	cdef int jx,jy
	
	cdef double r0,phi0
	cdef binning_theta = np.pi/20
	
	r0 = sqrt(x0**2+y0**2)
	if (r0==0.):
		phi0 = 0.0
	else:
		phi0 = acos(x0/r0)
		if (y0<0): phi0 = -phi0
	# fixing 2*n+abs(l) <= l_max
	
	
	ms = np.arange(-m_max,m_max+1,1,dtype=DTYPE2)
	creal = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
	cimag = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
	Mreal = np.zeros((n_max+1,2*l_max+1,2*m_max+1,l_max+1,l_max+1), dtype=DTYPE)
	Mimag = np.zeros((n_max+1,2*l_max+1,2*m_max+1,l_max+1,l_max+1), dtype=DTYPE)
	Mtotreal = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
	Mtotimag = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
 	
	####################################################################
	# scheme for M(t=0) and c
	for m in range(0,2*m_max+1):
		creal[0,0,m] =  psi(r0, phi0, theta0, d, 0, 0, ms[m])[0]
		cimag[0,0,m] = -psi(r0, phi0, theta0, d, 0, 0, ms[m])[1]
		Mreal[0,0,m,0,0] = creal[0,0,m]
		Mimag[0,0,m,0,0] = cimag[0,0,m]
		
		for l in range(1,l_max+1):
			for q in range(1,l+1):
				Mreal[0,l,m,q,q] = a* sqrt(1.0*l) * Mreal[0,l-1,m,q-1,q-1] / (landa(D_rot,mu,k,0,l,ms[m]) - landa(D_rot,mu,k,0,l-q,ms[m]))
				Mimag[0,l,m,q,q] = a* sqrt(1.0*l) * Mimag[0,l-1,m,q-1,q-1] / (landa(D_rot,mu,k,0,l,ms[m]) - landa(D_rot,mu,k,0,l-q,ms[m]))
				Mtotreal[0,l,m] += Mreal[0,l,m,q,q]
				Mtotimag[0,l,m] += Mimag[0,l,m,q,q]
				Mreal[0,-l,m,q,0] = a* sqrt(1.0*l) * Mreal[0,-l+1,m,q-1,0] / (landa(D_rot,mu,k,0,-l,ms[m]) - landa(D_rot,mu,k,0,-l+q,ms[m]))
				Mimag[0,-l,m,q,0] = a* sqrt(1.0*l) * Mimag[0,-l+1,m,q-1,0] / (landa(D_rot,mu,k,0,-l,ms[m]) - landa(D_rot,mu,k,0,-l+q,ms[m]))
				Mtotreal[0,-l,m] += Mreal[0,-l,m,q,0]
				Mtotimag[0,-l,m] += Mimag[0,-l,m,q,0]
			creal[0,l,m] =  psi(r0, phi0, theta0, d, 0, l, ms[m])[0] - Mtotreal[0,l,m]
			cimag[0,l,m] = -psi(r0, phi0, theta0, d, 0, l, ms[m])[1] - Mtotimag[0,l,m]
			Mreal[0,l,m,0,0] = creal[0,l,m]
			Mimag[0,l,m,0,0] = cimag[0,l,m]
			creal[0,-l,m] =  psi(r0, phi0, theta0, d, 0, -l, ms[m])[0] - Mtotreal[0,-l,m]
			cimag[0,-l,m] = -psi(r0, phi0, theta0, d, 0, -l, ms[m])[1] - Mtotimag[0,-l,m]
			Mreal[0,-l,m,0,0] = creal[0,-l,m]
			Mimag[0,-l,m,0,0] = cimag[0,-l,m]
				
		for n in range(1,n_max+1):
			for q in range(1,2*n+1):
				if (q<=n):
					Mreal[n,0,m,q,0] = -a* sqrt(1.0*n) * Mreal[n-1,1,m,q-1,0] / (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-q,q,ms[m]))
					Mimag[n,0,m,q,0] = -a* sqrt(1.0*n) * Mimag[n-1,1,m,q-1,0] / (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-q,q,ms[m]))
					Mtotreal[n,0,m] += Mreal[n,0,m,q,0]
					Mtotimag[n,0,m] += Mimag[n,0,m,q,0]
#~ 					if (n==1): print n,0,ms[m],q,0,Mreal[n,0,m,q,0]
				for p in range(max(1,q-n),min(n,q)+1):
					Mreal[n,0,m,q,p] = -a* sqrt(1.0*n) * (Mreal[n-1,1,m,q-1,p] + Mreal[n-1,-1,m,q-1,p-1])/ (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-(q+abs(q-2*p))/2,q-2*p,ms[m]))
					Mimag[n,0,m,q,p] = -a* sqrt(1.0*n) * (Mimag[n-1,1,m,q-1,p] + Mimag[n-1,-1,m,q-1,p-1])/ (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-(q+abs(q-2*p))/2,q-2*p,ms[m]))
					Mtotreal[n,0,m] += Mreal[n,0,m,q,p]
					Mtotimag[n,0,m] += Mimag[n,0,m,q,p]
#~ 					if (n==1): print n,0,ms[m],q,p,Mreal[n,0,m,q,p],Mreal[n-1,1,m,q-1,p],Mreal[n-1,-1,m,q-1,p-1]
			creal[n,0,m] =  psi(r0, phi0, theta0, d, n, 0, ms[m])[0] - Mtotreal[n,0,m]
			cimag[n,0,m] = -psi(r0, phi0, theta0, d, n, 0, ms[m])[1] - Mtotimag[n,0,m]
			Mreal[n,0,m,0,0] = creal[n,0,m]
			Mimag[n,0,m,0,0] = cimag[n,0,m]
			

			for l in range(1,l_max -2*n +1):
				for q in range(1,2*n+l+1):
					s = n-(q-l+abs(l+q))/2
					if (s>=0):
						Mreal[n,l,m,q,0] = -a* sqrt(1.0*n) * Mreal[n-1,l+1,m,q-1,0] / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q,ms[m]))
						Mimag[n,l,m,q,0] = -a* sqrt(1.0*n) * Mimag[n-1,l+1,m,q-1,0] / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q,ms[m]))
						Mtotreal[n,l,m] += Mreal[n,l,m,q,0]
						Mtotimag[n,l,m] += Mimag[n,l,m,q,0]
					for p in range(1,q+1):
						s = n-(q-l+abs(l+q-2*p))/2
						if (s>=0):
							Mreal[n,l,m,q,p] = a*(sqrt(1.0*(n+l))*Mreal[n,l-1,m,q-1,p-1] - sqrt(1.0*n)*Mreal[n-1,l+1,m,q-1,p]) / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q-2*p,ms[m]))
							Mimag[n,l,m,q,p] = a*(sqrt(1.0*(n+l))*Mimag[n,l-1,m,q-1,p-1] - sqrt(1.0*n)*Mimag[n-1,l+1,m,q-1,p]) / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q-2*p,ms[m]))
							Mtotreal[n,l,m] += Mreal[n,l,m,q,p]
							Mtotimag[n,l,m] += Mimag[n,l,m,q,p]
					
					s= n-(q-l+abs(-l+q))/2
					if (s>=0):
						Mreal[n,-l,m,q,0] = a* sqrt(1.0*(n+l)) * Mreal[n,-l+1,m,q-1,0] / (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q,ms[m]))
						Mimag[n,-l,m,q,0] = a* sqrt(1.0*(n+l)) * Mimag[n,-l+1,m,q-1,0] / (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q,ms[m]))
						Mtotreal[n,-l,m] += Mreal[n,-l,m,q,0]
						Mtotimag[n,-l,m] += Mimag[n,-l,m,q,0]
					for p in range(1,q+1):	
						s= n-(q-l+abs(-l+q-2*p))/2
						if (s>=0):
							Mreal[n,-l,m,q,p] = a* (sqrt(1.0*(n+l))*Mreal[n,-l+1,m,q-1,p] - sqrt(1.0*n)*Mreal[n-1,-l-1,m,q-1,p-1] )/ (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q-2*p,ms[m]))
							Mimag[n,-l,m,q,p] = a* (sqrt(1.0*(n+l))*Mimag[n,-l+1,m,q-1,p] - sqrt(1.0*n)*Mimag[n-1,-l-1,m,q-1,p-1] )/ (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q-2*p,ms[m]))
							Mtotreal[n,-l,m] += Mreal[n,-l,m,q,p]
							Mtotimag[n,-l,m] += Mimag[n,-l,m,q,p]
				creal[n,l,m] =  psi(r0, phi0, theta0, d, n, l, ms[m])[0] - Mtotreal[n,l,m]
				cimag[n,l,m] = -psi(r0, phi0, theta0, d, n, l, ms[m])[1] - Mtotimag[n,l,m]
				Mreal[n,l,m,0,0] = creal[n,l,m]
				Mimag[n,l,m,0,0] = cimag[n,l,m]
				creal[n,-l,m] =  psi(r0, phi0, theta0, d, n, -l, ms[m])[0] - Mtotreal[n,-l,m]
				cimag[n,-l,m] = -psi(r0, phi0, theta0, d, n, -l, ms[m])[1] - Mtotimag[n,-l,m]
				Mreal[n,-l,m,0,0] = creal[n,-l,m]
				Mimag[n,-l,m,0,0] = cimag[n,-l,m]
				
#~ 	print '---------------------------------'

	# scheme for M
	Mreal = np.zeros((n_max+1,2*l_max+1,2*m_max+1,l_max+1,l_max+1), dtype=DTYPE)
	Mimag = np.zeros((n_max+1,2*l_max+1,2*m_max+1,l_max+1,l_max+1), dtype=DTYPE)
	Mtotreal = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
	Mtotimag = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
	for m in range(0,2*m_max+1):
		Mreal[0,0,m,0,0] = creal[0,0,m] * exp(-landa(D_rot, mu, k, 0, 0, ms[m]) * t)
		Mimag[0,0,m,0,0] = cimag[0,0,m] * exp(-landa(D_rot, mu, k, 0, 0, ms[m]) * t)
		Mtotreal[0,0,m] = Mreal[0,0,m,0,0]
		Mtotimag[0,0,m] = Mimag[0,0,m,0,0]
#~ 		print 0,0,ms[m],0,0,Mreal[0,0,m,0,0]
		
		for l in range(1,l_max+1):
			Mreal[0,l,m,0,0] = creal[0,l,m] * exp(-landa(D_rot, mu, k, 0, l, ms[m]) * t)
			Mimag[0,l,m,0,0] = cimag[0,l,m] * exp(-landa(D_rot, mu, k, 0, l, ms[m]) * t)
			Mtotreal[0,l,m] = Mreal[0,l,m,0,0]
			Mtotimag[0,l,m] = Mimag[0,l,m,0,0]
			Mreal[0,-l,m,0,0] = creal[0,-l,m] * exp(-landa(D_rot, mu, k, 0, -l, ms[m]) * t)
			Mimag[0,-l,m,0,0] = cimag[0,-l,m] * exp(-landa(D_rot, mu, k, 0, -l, ms[m]) * t)
			Mtotreal[0,-l,m] = Mreal[0,-l,m,0,0]
			Mtotimag[0,-l,m] = Mimag[0,-l,m,0,0]
#~ 			print 0,l,ms[m],0,0,Mreal[0,l,m,0,0]
			for q in range(1,l+1):
				Mreal[0,l,m,q,q] = a* sqrt(1.0*l) * Mreal[0,l-1,m,q-1,q-1] / (landa(D_rot,mu,k,0,l,ms[m]) - landa(D_rot,mu,k,0,l-q,ms[m]))
				Mimag[0,l,m,q,q] = a* sqrt(1.0*l) * Mimag[0,l-1,m,q-1,q-1] / (landa(D_rot,mu,k,0,l,ms[m]) - landa(D_rot,mu,k,0,l-q,ms[m]))
				Mtotreal[0,l,m] += Mreal[0,l,m,q,q]
				Mtotimag[0,l,m] += Mimag[0,l,m,q,q]
				Mreal[0,-l,m,q,0] = a* sqrt(1.0*l) * Mreal[0,-l+1,m,q-1,0] / (landa(D_rot,mu,k,0,-l,ms[m]) - landa(D_rot,mu,k,0,-l+q,ms[m]))
				Mimag[0,-l,m,q,0] = a* sqrt(1.0*l) * Mimag[0,-l+1,m,q-1,0] / (landa(D_rot,mu,k,0,-l,ms[m]) - landa(D_rot,mu,k,0,-l+q,ms[m]))
				Mtotreal[0,-l,m] += Mreal[0,-l,m,q,0]
				Mtotimag[0,-l,m] += Mimag[0,-l,m,q,0]
#~ 				print 0,l,ms[m],q,q,Mreal[0,l,m,q,q]
				
		for n in range(1,n_max+1):
			Mreal[n,0,m,0,0] = creal[n,0,m] * exp(-landa(D_rot, mu, k, n, 0, ms[m]) * t)
			Mimag[n,0,m,0,0] = cimag[n,0,m] * exp(-landa(D_rot, mu, k, n, 0, ms[m]) * t)
			Mtotreal[n,0,m] = Mreal[n,0,m,0,0]
			Mtotimag[n,0,m] = Mimag[n,0,m,0,0]
#~ 			print n,0,ms[m],0,0,Mreal[n,0,m,0,0]
			for q in range(1,2*n+1):
				if (q<=n):
					Mreal[n,0,m,q,0] = -a* sqrt(1.0*n) * Mreal[n-1,1,m,q-1,0]/ (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-q,q,ms[m]))
					Mimag[n,0,m,q,0] = -a* sqrt(1.0*n) * Mimag[n-1,1,m,q-1,0]/ (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-q,q,ms[m]))
					Mtotreal[n,0,m] += Mreal[n,0,m,q,0]
					Mtotimag[n,0,m] += Mimag[n,0,m,q,0]
#~ 					print n,0,ms[m],q,0,Mreal[n,0,m,q,0]
				for p in range(max(1,q-n),min(n,q)+1):
					Mreal[n,0,m,q,p] = -a* sqrt(1.0*n) * (Mreal[n-1,1,m,q-1,p] + Mreal[n-1,-1,m,q-1,p-1])/ (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-(q+abs(q-2*p))/2,q-2*p,ms[m]))
					Mimag[n,0,m,q,p] = -a* sqrt(1.0*n) * (Mimag[n-1,1,m,q-1,p] + Mimag[n-1,-1,m,q-1,p-1])/ (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-(q+abs(q-2*p))/2,q-2*p,ms[m]))
					Mtotreal[n,0,m] += Mreal[n,0,m,q,p]
					Mtotimag[n,0,m] += Mimag[n,0,m,q,p]
#~ 					print n,0,ms[m],q,p,Mreal[n,0,m,q,p]

			for l in range(1,l_max -2*n +1):
				Mreal[n,l,m,0,0] = creal[n,l,m] * exp(-landa(D_rot, mu, k, n, l, ms[m]) * t)
				Mimag[n,l,m,0,0] = cimag[n,l,m] * exp(-landa(D_rot, mu, k, n, l, ms[m]) * t)
				Mtotreal[n,l,m] = Mreal[n,l,m,0,0]
				Mtotimag[n,l,m] = Mimag[n,l,m,0,0]
				Mreal[n,-l,m,0,0] = creal[n,-l,m] * exp(-landa(D_rot, mu, k, n, -l, ms[m]) * t)
				Mimag[n,-l,m,0,0] = cimag[n,-l,m] * exp(-landa(D_rot, mu, k, n, -l, ms[m]) * t)
				Mtotreal[n,-l,m] = Mreal[n,-l,m,0,0]
				Mtotimag[n,-l,m] = Mimag[n,-l,m,0,0]
#~ 				print n,l,ms[m],0,0,Mreal[n,l,m,0,0]
				
				for q in range(1,2*n+l+1):
					s = n-(q-l+abs(l+q))/2
					if (s>=0):
						Mreal[n,l,m,q,0] = -a* sqrt(1.0*n) * Mreal[n-1,l+1,m,q-1,0] / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q,ms[m]))
						Mimag[n,l,m,q,0] = -a* sqrt(1.0*n) * Mimag[n-1,l+1,m,q-1,0] / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q,ms[m]))
						Mtotreal[n,l,m] += Mreal[n,l,m,q,0]
						Mtotimag[n,l,m] += Mimag[n,l,m,q,0]
#~ 						print n,l,ms[m],q,0,Mreal[n,l,m,q,0]
					for p in range(1,q+1):
						s = n-(q-l+abs(l+q-2*p))/2
						if (s>=0):
							Mreal[n,l,m,q,p] = a*(sqrt(1.0*(n+l))*Mreal[n,l-1,m,q-1,p-1] - sqrt(1.0*n)*Mreal[n-1,l+1,m,q-1,p]) / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q-2*p,ms[m]))
							Mimag[n,l,m,q,p] = a*(sqrt(1.0*(n+l))*Mimag[n,l-1,m,q-1,p-1] - sqrt(1.0*n)*Mimag[n-1,l+1,m,q-1,p]) / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q-2*p,ms[m]))
							Mtotreal[n,l,m] += Mreal[n,l,m,q,p]
							Mtotimag[n,l,m] += Mimag[n,l,m,q,p]
#~ 							print n,l,ms[m],q,p,Mreal[n,l,m,q,p]
					
					s= n-(q-l+abs(-l+q))/2
					if (s>=0):
						Mreal[n,-l,m,q,0] = a* sqrt(1.0*(n+l)) * Mreal[n,-l+1,m,q-1,0] / (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q,ms[m]))
						Mimag[n,-l,m,q,0] = a* sqrt(1.0*(n+l)) * Mimag[n,-l+1,m,q-1,0] / (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q,ms[m]))
						Mtotreal[n,-l,m] += Mreal[n,-l,m,q,0]
						Mtotimag[n,-l,m] += Mimag[n,-l,m,q,0]
					for p in range(1,q+1):	
						s= n-(q-l+abs(-l+q-2*p))/2
						if (s>=0):
							Mreal[n,-l,m,q,p] = a* (sqrt(1.0*(n+l))*Mreal[n,-l+1,m,q-1,p] - sqrt(1.0*n)*Mreal[n-1,-l-1,m,q-1,p-1] )/ (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q-2*p,ms[m]))
							Mimag[n,-l,m,q,p] = a* (sqrt(1.0*(n+l))*Mimag[n,-l+1,m,q-1,p] - sqrt(1.0*n)*Mimag[n-1,-l-1,m,q-1,p-1] )/ (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q-2*p,ms[m]))
							Mtotreal[n,-l,m] += Mreal[n,-l,m,q,p]
							Mtotimag[n,-l,m] += Mimag[n,-l,m,q,p]
							
	for m in range(0,2*m_max+1):
		for n in range(0,n_max+1):
			if (ms[m]==0): print n,0,ms[m],Mtotreal[n,0,m]				
			for l in range(1,l_max -2*n +1):	
				if (ms[m]==0): print n,l,ms[m],Mtotreal[n,l,m]
				
#~ 	for m in range(0,2*m_max+1):
#~ 		for n in range(0,n_max+1):
#~ 			for q in range(0,2*n+1):
#~ 				for p in range(0,q+1):
#~ 					if (ms[m]==0): print n,0,ms[m],q,p,Mreal[n,0,m,q,p]
								
#~ 			for l in range(1,l_max -2*n +1):	
#~ 				for q in range(0,2*n+l+1):
#~ 					for p in range(0,q+1):
#~ 						if (ms[m]==0): print n,l,ms[m],q,p,Mreal[n,l,m,q,p]
	
	xs = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning)
	ys = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning)
#~ 	thetas = np.arange(0,2*np.pi,binning_theta)
	data = np.zeros((Nbins,Nbins), dtype=DTYPE)
	
	norma = 0.0
	for jx in range(0,len(xs)):
		for jy in range(0,len(ys)):
			x=xs[jx]
			y=ys[jy]
			r = sqrt(x**2+y**2)
			if (r==0.):
				phi=0.0
			else:
				phi = acos(x/r)
				if (y<0): phi = -phi

			sumareal = 0.0
			for n in range(0,n_max+1):
				sumareal +=  Mtotreal[n,0,m_max] * psi_integrated_theta (r, phi, d, n, 0)[0] 
#~ 				sumareal += -Mtotimag[n,0,m_max] * psi_integrated_theta (r, phi, d, n, 0)[1]
				for l in range(1,l_max -2*n +1):
					m = l + m_max
					sumareal +=  2*Mtotreal[n,l,m] * psi_integrated_theta (r, phi, d, n, l)[0] 
			#		sumareal +=  Mtotreal[n,-l,-l+m_max] * psi_integrated_theta (r, phi, d, n, -l)[0]
					sumareal += -2*Mtotimag[n,l,m] * psi_integrated_theta (r, phi, d, n, l)[1] 
#~ 					sumareal +=  -Mtotimag[n,-l,-l+m_max] * psi_integrated_theta (r, phi, d, n, -l)[1]
#~ 					if (jx==0 and jy==0): print n,l,ms[m],ms[-l+m_max],Mtotreal[n,l,m],Mtotreal[n,-l,-l+m_max],Mtotimag[n,l,m],Mtotimag[n,-l,-l+m_max]

			sumareal = sumareal* peq (r, d)
			data[jx,jy] = sumareal
#~ 			if (sumareal<0.0): data[jx,jy] = 0.0
			norma += sumareal*binning**2
	print 'norma = ', norma

	return data

def spatial_distribution_numerics_v2(double D_rot, double D, double v, double k, double mu, double x0, double y0, double theta0, double t, double binning, int Nbins, int qmax):		
	# see main.py for a definition of various parameters
	cdef int n,l,m
	cdef double x,y,theta
	cdef double r,phi
	cdef double d = mu*k/(2*D)
	cdef double a = v*sqrt(d)
	cdef int n_max = qmax/2
	cdef int l_max=2*n_max
	cdef int m_max = l_max
#~ 	cdef int m_max = 0
	# m_max should be equal to l_max
	cdef int q,p,s,l2
	cdef int jx,jy
	
	cdef double r0,phi0
	cdef binning_theta = np.pi/20
	
	r0 = sqrt(x0**2+y0**2)
	if (r0==0.):
		phi0 = 0.0
	else:
		phi0 = acos(x0/r0)
		if (y0<0): phi0 = -phi0
	# fixing 2*n+abs(l) <= l_max
	
	
	ms = np.arange(-m_max,m_max+1,1,dtype=DTYPE2)
	creal = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
	cimag = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
	Mreal0 = np.zeros((n_max+1,2*l_max+1,2*m_max+1,l_max+1,l_max+1), dtype=DTYPE)
	Mimag0 = np.zeros((n_max+1,2*l_max+1,2*m_max+1,l_max+1,l_max+1), dtype=DTYPE)
	Mtotreal = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
	Mtotimag = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
 	
	####################################################################
	# scheme for M(t=0) and c
	for m in range(0,2*m_max+1):
		creal[0,0,m] =  psi(r0, phi0, theta0, d, 0, 0, ms[m])[0]
		cimag[0,0,m] = -psi(r0, phi0, theta0, d, 0, 0, ms[m])[1]
		Mreal0[0,0,m,0,0] = creal[0,0,m]
		Mimag0[0,0,m,0,0] = cimag[0,0,m]
		
		for l in range(1,l_max+1):
			for q in range(1,l+1):
				Mreal0[0,l,m,q,q] = a* sqrt(1.0*l) * Mreal0[0,l-1,m,q-1,q-1] / (landa(D_rot,mu,k,0,l,ms[m]) - landa(D_rot,mu,k,0,l-q,ms[m]))
				Mimag0[0,l,m,q,q] = a* sqrt(1.0*l) * Mimag0[0,l-1,m,q-1,q-1] / (landa(D_rot,mu,k,0,l,ms[m]) - landa(D_rot,mu,k,0,l-q,ms[m]))
				Mtotreal[0,l,m] += Mreal0[0,l,m,q,q]
				Mtotimag[0,l,m] += Mimag0[0,l,m,q,q]
				Mreal0[0,-l,m,q,0] = a* sqrt(1.0*l) * Mreal0[0,-l+1,m,q-1,0] / (landa(D_rot,mu,k,0,-l,ms[m]) - landa(D_rot,mu,k,0,-l+q,ms[m]))
				Mimag0[0,-l,m,q,0] = a* sqrt(1.0*l) * Mimag0[0,-l+1,m,q-1,0] / (landa(D_rot,mu,k,0,-l,ms[m]) - landa(D_rot,mu,k,0,-l+q,ms[m]))
				Mtotreal[0,-l,m] += Mreal0[0,-l,m,q,0]
				Mtotimag[0,-l,m] += Mimag0[0,-l,m,q,0]
			creal[0,l,m] =  psi(r0, phi0, theta0, d, 0, l, ms[m])[0] - Mtotreal[0,l,m]
			cimag[0,l,m] = -psi(r0, phi0, theta0, d, 0, l, ms[m])[1] - Mtotimag[0,l,m]
			Mreal0[0,l,m,0,0] = creal[0,l,m]
			Mimag0[0,l,m,0,0] = cimag[0,l,m]
			creal[0,-l,m] =  psi(r0, phi0, theta0, d, 0, -l, ms[m])[0] - Mtotreal[0,-l,m]
			cimag[0,-l,m] = -psi(r0, phi0, theta0, d, 0, -l, ms[m])[1] - Mtotimag[0,-l,m]
			Mreal0[0,-l,m,0,0] = creal[0,-l,m]
			Mimag0[0,-l,m,0,0] = cimag[0,-l,m]
				
		for n in range(1,n_max+1):
			for q in range(1,2*n+1):
				if (q<=n):
					Mreal0[n,0,m,q,0] = -a* sqrt(1.0*n) * Mreal0[n-1,1,m,q-1,0] / (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-q,q,ms[m]))
					Mimag0[n,0,m,q,0] = -a* sqrt(1.0*n) * Mimag0[n-1,1,m,q-1,0] / (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-q,q,ms[m]))
					Mtotreal[n,0,m] += Mreal0[n,0,m,q,0]
					Mtotimag[n,0,m] += Mimag0[n,0,m,q,0]
#~ 					if (n==1): print n,0,ms[m],q,0,Mreal[n,0,m,q,0]
				for p in range(max(1,q-n),min(n,q)+1):
					Mreal0[n,0,m,q,p] = -a* sqrt(1.0*n) * (Mreal0[n-1,1,m,q-1,p] + Mreal0[n-1,-1,m,q-1,p-1])/ (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-(q+abs(q-2*p))/2,q-2*p,ms[m]))
					Mimag0[n,0,m,q,p] = -a* sqrt(1.0*n) * (Mimag0[n-1,1,m,q-1,p] + Mimag0[n-1,-1,m,q-1,p-1])/ (landa(D_rot,mu,k,n,0,ms[m]) - landa(D_rot,mu,k,n-(q+abs(q-2*p))/2,q-2*p,ms[m]))
					Mtotreal[n,0,m] += Mreal0[n,0,m,q,p]
					Mtotimag[n,0,m] += Mimag0[n,0,m,q,p]
#~ 					if (n==1): print n,0,ms[m],q,p,Mreal[n,0,m,q,p],Mreal[n-1,1,m,q-1,p],Mreal[n-1,-1,m,q-1,p-1]
			creal[n,0,m] =  psi(r0, phi0, theta0, d, n, 0, ms[m])[0] - Mtotreal[n,0,m]
			cimag[n,0,m] = -psi(r0, phi0, theta0, d, n, 0, ms[m])[1] - Mtotimag[n,0,m]
			Mreal0[n,0,m,0,0] = creal[n,0,m]
			Mimag0[n,0,m,0,0] = cimag[n,0,m]
			

			for l in range(1,l_max -2*n +1):
				for q in range(1,2*n+l+1):
					s = n-(q-l+abs(l+q))/2
					if (s>=0):
						Mreal0[n,l,m,q,0] = -a* sqrt(1.0*n) * Mreal0[n-1,l+1,m,q-1,0] / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q,ms[m]))
						Mimag0[n,l,m,q,0] = -a* sqrt(1.0*n) * Mimag0[n-1,l+1,m,q-1,0] / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q,ms[m]))
						Mtotreal[n,l,m] += Mreal0[n,l,m,q,0]
						Mtotimag[n,l,m] += Mimag0[n,l,m,q,0]
					for p in range(1,q+1):
						s = n-(q-l+abs(l+q-2*p))/2
						if (s>=0):
							Mreal0[n,l,m,q,p] = a*(sqrt(1.0*(n+l))*Mreal0[n,l-1,m,q-1,p-1] - sqrt(1.0*n)*Mreal0[n-1,l+1,m,q-1,p]) / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q-2*p,ms[m]))
							Mimag0[n,l,m,q,p] = a*(sqrt(1.0*(n+l))*Mimag0[n,l-1,m,q-1,p-1] - sqrt(1.0*n)*Mimag0[n-1,l+1,m,q-1,p]) / (landa(D_rot,mu,k,n,l,ms[m]) - landa(D_rot,mu,k,s,l+q-2*p,ms[m]))
							Mtotreal[n,l,m] += Mreal0[n,l,m,q,p]
							Mtotimag[n,l,m] += Mimag0[n,l,m,q,p]
					
					s= n-(q-l+abs(-l+q))/2
					if (s>=0):
						Mreal0[n,-l,m,q,0] = a* sqrt(1.0*(n+l)) * Mreal0[n,-l+1,m,q-1,0] / (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q,ms[m]))
						Mimag0[n,-l,m,q,0] = a* sqrt(1.0*(n+l)) * Mimag0[n,-l+1,m,q-1,0] / (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q,ms[m]))
						Mtotreal[n,-l,m] += Mreal0[n,-l,m,q,0]
						Mtotimag[n,-l,m] += Mimag0[n,-l,m,q,0]
					for p in range(1,q+1):	
						s= n-(q-l+abs(-l+q-2*p))/2
						if (s>=0):
							Mreal0[n,-l,m,q,p] = a* (sqrt(1.0*(n+l))*Mreal0[n,-l+1,m,q-1,p] - sqrt(1.0*n)*Mreal0[n-1,-l-1,m,q-1,p-1] )/ (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q-2*p,ms[m]))
							Mimag0[n,-l,m,q,p] = a* (sqrt(1.0*(n+l))*Mimag0[n,-l+1,m,q-1,p] - sqrt(1.0*n)*Mimag0[n-1,-l-1,m,q-1,p-1] )/ (landa(D_rot,mu,k,n,-l,ms[m]) - landa(D_rot,mu,k,s,-l+q-2*p,ms[m]))
							Mtotreal[n,-l,m] += Mreal0[n,-l,m,q,p]
							Mtotimag[n,-l,m] += Mimag0[n,-l,m,q,p]
				creal[n,l,m] =  psi(r0, phi0, theta0, d, n, l, ms[m])[0] - Mtotreal[n,l,m]
				cimag[n,l,m] = -psi(r0, phi0, theta0, d, n, l, ms[m])[1] - Mtotimag[n,l,m]
				Mreal0[n,l,m,0,0] = creal[n,l,m]
				Mimag0[n,l,m,0,0] = cimag[n,l,m]
				creal[n,-l,m] =  psi(r0, phi0, theta0, d, n, -l, ms[m])[0] - Mtotreal[n,-l,m]
				cimag[n,-l,m] = -psi(r0, phi0, theta0, d, n, -l, ms[m])[1] - Mtotimag[n,-l,m]
				Mreal0[n,-l,m,0,0] = creal[n,-l,m]
				Mimag0[n,-l,m,0,0] = cimag[n,-l,m]
				
#~ 	print '---------------------------------'

	# scheme for M
	Mtotreal = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
	Mtotimag = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)		
	for m in range(0,2*m_max+1):
		for n in range(0,n_max+1):
			Mtotreal[n,0,m]+= psi(r0, phi0, theta0, d, n, 0, ms[m])[0] * exp(-landa(D_rot, mu, k, n, 0, ms[m]) * t)
			Mtotimag[n,0,m]-= psi(r0, phi0, theta0, d, n, 0, ms[m])[1] * exp(-landa(D_rot, mu, k, n, 0, ms[m]) * t)
			for q in range(1,2*n+1):
				for p in range(0,q+1):
					s = n-(q+abs(q-2*p))/2
					l2 = q-2*p
					if (s>=0):
						Mtotreal[n,0,m]+=Mreal0[n,0,m,q,p]*(exp(-landa(D_rot, mu, k, s, l2, ms[m]) * t)-exp(-landa(D_rot, mu, k, n, 0, ms[m]) * t))
						Mtotimag[n,0,m]+=Mimag0[n,0,m,q,p]*(exp(-landa(D_rot, mu, k, s, l2, ms[m]) * t)-exp(-landa(D_rot, mu, k, n, 0, ms[m]) * t))
#~ 						if (ms[m]==0): print n,0,ms[m],q,p,Mreal[n,0,m,q,p],Mreal0[n,0,m,q,p],Mreal[n,0,m,q,p]-Mreal0[n,0,m,q,p]*exp(-landa(D_rot, mu, k, s, l2, ms[m]) * t)
								
			for l in range(1,l_max -2*n +1):
				Mtotreal[n,l,m]+= psi(r0, phi0, theta0, d, n, l, ms[m])[0] * exp(-landa(D_rot, mu, k, n, l, ms[m]) * t)
				Mtotimag[n,l,m]-= psi(r0, phi0, theta0, d, n, l, ms[m])[1] * exp(-landa(D_rot, mu, k, n, l, ms[m]) * t)
				for q in range(1,2*n+l+1):
					for p in range(0,q+1):
						s = n-(q-l+abs(l+q-2*p))/2
						l2 = l+q-2*p
						if (s>=0):
							Mtotreal[n,l,m]+=Mreal0[n,l,m,q,p]*(exp(-landa(D_rot, mu, k, s, l2, ms[m]) * t)-exp(-landa(D_rot, mu, k, n, l, ms[m]) * t))
							Mtotimag[n,l,m]+=Mimag0[n,l,m,q,p]*(exp(-landa(D_rot, mu, k, s, l2, ms[m]) * t)-exp(-landa(D_rot, mu, k, n, l, ms[m]) * t))
#~ 							if (ms[m]==0): print n,l,ms[m],q,p,Mreal[n,l,m,q,p],Mreal0[n,l,m,q,p],Mreal[n,l,m,q,p]-Mreal0[n,l,m,q,p]*exp(-landa(D_rot, mu, k, s, l2, ms[m]) * t)
	
	xs = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning)
	ys = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning)
#~ 	thetas = np.arange(0,2*np.pi,binning_theta)
	data = np.zeros((Nbins,Nbins), dtype=DTYPE)
	
	norma = 0.0
	for jx in range(0,len(xs)):
		for jy in range(0,len(ys)):
			x=xs[jx]
			y=ys[jy]
			r = sqrt(x**2+y**2)
			if (r==0.):
				phi=0.0
			else:
				phi = acos(x/r)
				if (y<0): phi = -phi

			sumareal = 0.0
			for n in range(0,n_max+1):
				sumareal +=  Mtotreal[n,0,m_max] * psi_integrated_theta (r, phi, d, n, 0)[0] 
#~ 				sumareal += -Mtotimag[n,0,m_max] * psi_integrated_theta (r, phi, d, n, 0)[1]
				for l in range(1,l_max -2*n +1):
					m = l + m_max
					sumareal +=  2*Mtotreal[n,l,m] * psi_integrated_theta (r, phi, d, n, l)[0] 
			#		sumareal +=  Mtotreal[n,-l,-l+m_max] * psi_integrated_theta (r, phi, d, n, -l)[0]
					sumareal += -2*Mtotimag[n,l,m] * psi_integrated_theta (r, phi, d, n, l)[1] 
#~ 					sumareal +=  -Mtotimag[n,-l,-l+m_max] * psi_integrated_theta (r, phi, d, n, -l)[1]
#~ 					if (jx==0 and jy==0): print n,l,ms[m],ms[-l+m_max],Mtotreal[n,l,m],Mtotreal[n,-l,-l+m_max],Mtotimag[n,l,m],Mtotimag[n,-l,-l+m_max]

			sumareal = sumareal* peq (r, d)
			data[jx,jy] = sumareal
#~ 			if (sumareal<0.0): data[jx,jy] = 0.0
			norma += sumareal*binning**2
	print 'norma = ', norma

	return data

def spatial_distribution_numerics(double D_rot, double D, double v, double k, double mu, double x0, double y0, double theta0, double t, double binning, int Nbins, int qmax):		
	# see main.py for a definition of various parameters
	cdef int n,l,m
	cdef double x,y,theta
	cdef double r,phi
	cdef double d = mu*k/(2*D)
	cdef double a = v*sqrt(d)
	cdef int n_max = qmax/2
	cdef int l_max=2*n_max
	cdef int m_max = l_max
#~ 	cdef int m_max = 0
	# m_max should be equal to l_max
	cdef int q,p,s,l2,lindex,q2,p2,n3,l3
	cdef double fp,fm
	cdef int jx,jy
	cdef double factor
	
	cdef double r0,phi0
	cdef binning_theta = np.pi/20
	
	r0 = sqrt(x0**2+y0**2)
	if (r0==0.):
		phi0 = 0.0
	else:
		phi0 = acos(x0/r0)
		if (y0<0): phi0 = -phi0
	# fixing 2*n+abs(l) <= l_max
	
	
	ms = np.arange(-m_max,m_max+1,1,dtype=DTYPE2)
	Creal = np.zeros((n_max+1,2*l_max+1,2*m_max+1,l_max+1,l_max+1), dtype=DTYPE)
 	
	####################################################################
	# scheme for Creal
	for m in range(0,2*m_max+1):
		n = 0
		l = 0
		Creal[n,l,m,0,0] = exp(-landa(D_rot, mu, k, n, l, ms[m]) * t)
		
		for lindex in range(1,l_max+1):
			l=lindex
			Creal[n,l,m,0,0] = exp(-landa(D_rot, mu, k, n, l, ms[m]) * t)
			l=-lindex
			Creal[n,l,m,0,0] = exp(-landa(D_rot, mu, k, n, l, ms[m]) * t)
			for q in range(1,lindex+1):
				l=lindex
				p=q
				l2 = l-q
				s=0
				fp = sqrt(1.*(l2+1))
				if (landa(D_rot, mu, k, n, l, ms[m]) != landa(D_rot, mu, k, s, l2, ms[m])):
					Creal[n,l,m,q,p] = (sqrt(1.*(n+l))*Creal[n,l-1,m,q-1,p-1] - fp*Creal[n,l,m,q-1,p-1])/(landa(D_rot, mu, k, n, l, ms[m])-landa(D_rot, mu, k, s, l2, ms[m]))
				else:
					Dreal = np.zeros((l_max+1,l_max+2), dtype=DTYPE)
					Dreal[1,1] = t * fp * Creal[s,l2,m,0,0]
					for q2 in range(2,q+1):
						p2=q2
						l3=l2-q2+2*p2
						n3=s+(q2-abs(l3)+abs(l2))/2
						if (l3>0):
							Dreal[q2,p2]  =  sqrt(1.*(n3+l3))*(Creal[n3,l3-1,m,q2-1,p2-1]-Dreal[q2-1,p2-1])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3, l3-1, ms[m]))
						else:
							Dreal[q2,p2]  = 0.0
					Creal[n,l,m,q,p] = Dreal[q,p]
				
				
				l=-lindex
				p=0	
				l2 = l+q
				s=0
				fm = sqrt(1.*(-l2+1))
				if (landa(D_rot, mu, k, n, l, ms[m]) != landa(D_rot, mu, k, s, l2, ms[m])):
					Creal[n,l,m,q,p] = (+sqrt(1.*(n-l))*Creal[n,l+1,m,q-1,p] - fm*Creal[n,l,m,q-1,p] )/(landa(D_rot, mu, k, n, l, ms[m])-landa(D_rot, mu, k, s, l2, ms[m]))
				else:
					Dreal = np.zeros((l_max+1,l_max+2), dtype=DTYPE)
					Dreal[1,0] = t * fm * Creal[s,l2,m,0,0]
					for q2 in range(2,q+1):
						p2=0
						l3=l2-q2+2*p2
						n3=s+(q2-abs(l3)+abs(l2))/2
						if (l3<0):
							Dreal[q2,p2] =  sqrt(1.*(n3-l3))*(Creal[n3,l3+1,m,q2-1,p2]-Dreal[q2-1,p2])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3, l3+1, ms[m]))
						else:
							Dreal[q2,p2] = 0.0
					Creal[n,l,m,q,p] = Dreal[q,p]
		
		for n in range(1,n_max+1):
			l = 0
			Creal[n,l,m,0,0] = exp(-landa(D_rot, mu, k, n, l, ms[m]) * t)
			for q in range(1,2*n+1):
				for p in range(0,q+1):
					l2 = q-2*p
					s=n-(q+abs(l2))/2
					if (s>=0):
						if (l2>=0):
							fp = sqrt(1.*(s+l2+1))
						else:
							fp = -sqrt(1.*(s+1))
						if (l2>0):
							fm = -sqrt(1.*(s+1))
						else:
							fm = sqrt(1.*(s-l2+1))
						if (landa(D_rot, mu, k, n, l, ms[m]) != landa(D_rot, mu, k, s, l2, ms[m])):
							Creal[n,l,m,q,p] = (-sqrt(1.*n)*Creal[n-1,l-1,m,q-1,p-1] - fp*Creal[n,l,m,q-1,p-1] -sqrt(1.*n)*Creal[n-1,l+1,m,q-1,p] - fm*Creal[n,l,m,q-1,p] )/(landa(D_rot, mu, k, n, l, ms[m])-landa(D_rot, mu, k, s, l2, ms[m]))
						else:
							if (ms[m]==0): print n,l,ms[m]
							if (ms[m]==0): print s,l2,q,p
							Dreal = np.zeros((l_max+1,l_max+2), dtype=DTYPE)
							Dreal[1,1] = t * fp * Creal[s,l2,m,0,0]
							Dreal[1,0] = t * fm * Creal[s,l2,m,0,0]
							for q2 in range(2,q+1):
								for p2 in range(0,min(q2+1,p+1)):
									l3=l2-q2+2*p2
									n3=s+(q2-abs(l3)+abs(l2))/2
									if (l3>0):
										Dreal[q2,p2]  =  sqrt(1.*(n3+l3))*(Creal[n3,l3-1,m,q2-1,p2-1]-Dreal[q2-1,p2-1])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3, l3-1, ms[m]))
										Dreal[q2,p2] += -sqrt(1.*n3)*     (Creal[n3-1,l3+1,m,q2-1,p2]-Dreal[q2-1,p2])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3+1, ms[m]))
									elif (l3==0):
										Dreal[q2,p2]  = -sqrt(1.*n3)*(Creal[n3-1,l3-1,m,q2-1,p2-1]-Dreal[q2-1,p2-1])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3-1, ms[m]))
										Dreal[q2,p2] += -sqrt(1.*n3)*(Creal[n3-1,l3+1,m,q2-1,p2]-Dreal[q2-1,p2])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3+1, ms[m]))
									else:
										Dreal[q2,p2]  = -sqrt(1.*n3)*(Creal[n3-1,l3-1,m,q2-1,p2-1]-Dreal[q2-1,p2-1])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3-1, ms[m]))
										Dreal[q2,p2] +=  sqrt(1.*(n3-l3))*(Creal[n3,l3+1,m,q2-1,p2]-Dreal[q2-1,p2])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3, l3+1, ms[m]))
									if (ms[m]==0): print 'here', q2,p2,n3,l3,Dreal[q2,p2]
							if (ms[m]==0): print Dreal[q,p]
							Creal[n,l,m,q,p] = Dreal[q,p]
			
			for lindex in range(1,l_max -2*n +1):
				l=lindex
				Creal[n,l,m,0,0] = exp(-landa(D_rot, mu, k, n, l, ms[m]) * t)
				l=-lindex
				Creal[n,l,m,0,0] = exp(-landa(D_rot, mu, k, n, l, ms[m]) * t)
				for q in range(1,2*n+lindex+1):
					for p in range(0,q+1):
						l=lindex
						l2 = l+q-2*p
						s=n-(q-l+abs(l2))/2
						if (s>=0):
							if (l2>=0):
								fp = sqrt(1.*(s+l2+1))
							else:
								fp = -sqrt(1.*(s+1))
							if (l2>0):
								fm = -sqrt(1.*(s+1))
							else:
								fm = sqrt(1.*(s-l2+1))
							if (landa(D_rot, mu, k, n, l, ms[m]) != landa(D_rot, mu, k, s, l2, ms[m])):
								Creal[n,l,m,q,p] = (sqrt(1.*(n+l))*Creal[n,l-1,m,q-1,p-1] - fp*Creal[n,l,m,q-1,p-1] -sqrt(1.*n)*Creal[n-1,l+1,m,q-1,p] - fm*Creal[n,l,m,q-1,p] )/(landa(D_rot, mu, k, n, l, ms[m])-landa(D_rot, mu, k, s, l2, ms[m]))
							else:
								Dreal = np.zeros((l_max+1,l_max+2), dtype=DTYPE)
								Dreal[1,1] = t * fp * Creal[s,l2,m,0,0]
								Dreal[1,0] = t * fm * Creal[s,l2,m,0,0]
								for q2 in range(2,q+1):
									for p2 in range(0,min(q2+1,p+1)):
										l3=l2-q2+2*p2
										n3=s+(q2-abs(l3)+abs(l2))/2
										if (l3>0):
											Dreal[q2,p2]  =  sqrt(1.*(n3+l3))*(Creal[n3,l3-1,m,q2-1,p2-1]-Dreal[q2-1,p2-1])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3, l3-1, ms[m]))
											Dreal[q2,p2] += -sqrt(1.*n3)*     (Creal[n3-1,l3+1,m,q2-1,p2]-Dreal[q2-1,p2])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3+1, ms[m]))
										elif (l3==0):
											Dreal[q2,p2]  = -sqrt(1.*n3)*(Creal[n3-1,l3-1,m,q2-1,p2-1]-Dreal[q2-1,p2-1])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3-1, ms[m]))
											Dreal[q2,p2] += -sqrt(1.*n3)*(Creal[n3-1,l3+1,m,q2-1,p2]-Dreal[q2-1,p2])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3+1, ms[m]))
										else:
											Dreal[q2,p2]  = -sqrt(1.*n3)*(Creal[n3-1,l3-1,m,q2-1,p2-1]-Dreal[q2-1,p2-1])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3-1, ms[m]))
											Dreal[q2,p2] +=  sqrt(1.*(n3-l3))*(Creal[n3,l3+1,m,q2-1,p2]-Dreal[q2-1,p2])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3, l3+1, ms[m]))
								Creal[n,l,m,q,p] = Dreal[q,p]
								
						l=-lindex
						l2 = l+q-2*p
						s=n-(q+l+abs(l2))/2
						if (s>=0):
							if (l2>=0):
								fp = sqrt(1.*(s+l2+1))
							else:
								fp = -sqrt(1.*(s+1))
							if (l2>0):
								fm = -sqrt(1.*(s+1))
							else:
								fm = sqrt(1.*(s-l2+1))
							if (landa(D_rot, mu, k, n, l, ms[m]) != landa(D_rot, mu, k, s, l2, ms[m])):
								Creal[n,l,m,q,p] = (-sqrt(1.*n)*Creal[n-1,l-1,m,q-1,p-1] - fp*Creal[n,l,m,q-1,p-1] +sqrt(1.*(n-l))*Creal[n,l+1,m,q-1,p] - fm*Creal[n,l,m,q-1,p] )/(landa(D_rot, mu, k, n, l, ms[m])-landa(D_rot, mu, k, s, l2, ms[m]))
							else:
								Dreal = np.zeros((l_max+1,l_max+2), dtype=DTYPE)
								Dreal[1,1] = t * fp * Creal[s,l2,m,0,0]
								Dreal[1,0] = t * fm * Creal[s,l2,m,0,0]
								for q2 in range(2,q+1):
									for p2 in range(0,min(q2+1,p+1)):
										l3=l2-q2+2*p2
										n3=s+(q2-abs(l3)+abs(l2))/2
										if (l3>0):
											Dreal[q2,p2]  =  sqrt(1.*(n3+l3))*(Creal[n3,l3-1,m,q2-1,p2-1]-Dreal[q2-1,p2-1])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3, l3-1, ms[m]))
											Dreal[q2,p2] += -sqrt(1.*n3)*     (Creal[n3-1,l3+1,m,q2-1,p2]-Dreal[q2-1,p2])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3+1, ms[m]))
										elif (l3==0):
											Dreal[q2,p2]  = -sqrt(1.*n3)*(Creal[n3-1,l3-1,m,q2-1,p2-1]-Dreal[q2-1,p2-1])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3-1, ms[m]))
											Dreal[q2,p2] += -sqrt(1.*n3)*(Creal[n3-1,l3+1,m,q2-1,p2]-Dreal[q2-1,p2])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3+1, ms[m]))
										else:
											Dreal[q2,p2]  = -sqrt(1.*n3)*(Creal[n3-1,l3-1,m,q2-1,p2-1]-Dreal[q2-1,p2-1])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3-1, l3-1, ms[m]))
											Dreal[q2,p2] +=  sqrt(1.*(n3-l3))*(Creal[n3,l3+1,m,q2-1,p2]-Dreal[q2-1,p2])/(landa(D_rot, mu, k, s, l2, ms[m])-landa(D_rot, mu, k, n3, l3+1, ms[m]))
								Creal[n,l,m,q,p] = Dreal[q,p]
							
				
#~ 	print '---------------------------------'

	# scheme for M
	Mtotreal = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)
	Mtotimag = np.zeros((n_max+1,2*l_max+1,2*m_max+1), dtype=DTYPE)		
	for m in range(0,2*m_max+1):
		for n in range(0,n_max+1):
			Mtotreal[n,0,m]+= psi(r0, phi0, theta0, d, n, 0, ms[m])[0] * Creal[n,0,m,0,0]
			Mtotimag[n,0,m]-= psi(r0, phi0, theta0, d, n, 0, ms[m])[1] * Creal[n,0,m,0,0]
			for q in range(1,2*n+1):
				for p in range(0,q+1):
					s = n-(q+abs(q-2*p))/2
					l2 = q-2*p
					if (s>=0):
						Mtotreal[n,0,m]+=a**q * Creal[n,0,m,q,p] * psi(r0, phi0, theta0, d, s, l2, ms[m])[0]
						Mtotimag[n,0,m]-=a**q * Creal[n,0,m,q,p] * psi(r0, phi0, theta0, d, s, l2, ms[m])[1]
			if (ms[m]==0): print n,0,ms[m],Mtotreal[n,0,m]
								
			for l in range(1,l_max -2*n +1):
				Mtotreal[n,l,m]+= psi(r0, phi0, theta0, d, n, l, ms[m])[0] * Creal[n,l,m,0,0]
				Mtotimag[n,l,m]-= psi(r0, phi0, theta0, d, n, l, ms[m])[1] * Creal[n,l,m,0,0]
				for q in range(1,2*n+l+1):
					for p in range(0,q+1):
						s = n-(q-l+abs(l+q-2*p))/2
						l2 = l+q-2*p
						if (s>=0):
							Mtotreal[n,l,m]+=a**q * Creal[n,l,m,q,p] * psi(r0, phi0, theta0, d, s, l2, ms[m])[0]
							Mtotimag[n,l,m]-=a**q * Creal[n,l,m,q,p] * psi(r0, phi0, theta0, d, s, l2, ms[m])[1]
				if (ms[m]==0): print n,l,ms[m],Mtotreal[n,l,m]
	
	xs = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning)
	ys = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning)
#~ 	thetas = np.arange(0,2*np.pi,binning_theta)
	data = np.zeros((Nbins,Nbins), dtype=DTYPE)
	
	norma = 0.0
	for jx in range(0,len(xs)):
		for jy in range(0,len(ys)):
			x=xs[jx]
			y=ys[jy]
			r = sqrt(x**2+y**2)
			if (r==0.):
				phi=0.0
			else:
				phi = acos(x/r)
				if (y<0): phi = -phi

			sumareal = 0.0
			for n in range(0,n_max+1):
				sumareal +=  Mtotreal[n,0,m_max] * psi_integrated_theta (r, phi, d, n, 0)[0] 
#~ 				sumareal += -Mtotimag[n,0,m_max] * psi_integrated_theta (r, phi, d, n, 0)[1]
				for l in range(1,l_max -2*n +1):
					m = l + m_max
					sumareal +=  2*Mtotreal[n,l,m] * psi_integrated_theta (r, phi, d, n, l)[0] 
			#		sumareal +=  Mtotreal[n,-l,-l+m_max] * psi_integrated_theta (r, phi, d, n, -l)[0]
					sumareal += -2*Mtotimag[n,l,m] * psi_integrated_theta (r, phi, d, n, l)[1] 
#~ 					sumareal +=  -Mtotimag[n,-l,-l+m_max] * psi_integrated_theta (r, phi, d, n, -l)[1]
#~ 					if (jx==0 and jy==0): print n,l,ms[m],ms[-l+m_max],Mtotreal[n,l,m],Mtotreal[n,-l,-l+m_max],Mtotimag[n,l,m],Mtotimag[n,-l,-l+m_max]

			sumareal = sumareal* peq (r, d)
			data[jx,jy] = sumareal
#~ 			if (sumareal<0.0): data[jx,jy] = 0.0
			norma += sumareal*binning**2
	print 'norma = ', norma

	return data