# nonlocaleigvalue.py
#  
# Compute dispersion relation from the transverse instability NLEP
# of an interior homoclinic stripe for several Ly lengths and x0.
#
#---------------
#
import os                         # For issuing commands to the OS.
import sys                        # For determining the Python version.
#
# Algebra and plotting libraries
import matplotlib
from matplotlib import rcParams
from pylab import *
from numpy import *
from scipy import *
from scipy import linalg
from scipy.linalg import eig
from scipy import sparse
#
rcParams['xtick.direction'] = 'out'
rcParams['ytick.direction'] = 'out'
#
# Clear command shell
clear = lambda: os.system('clear')
clear()
#
# Print some useful information
#
print 'Executing on', os.uname()
print 'Python version', sys.version
print 'matplotlib version', matplotlib.__version__
#
# Original parameter values
#
k1 = 0.008
b  = 0.008
c  = 0.1
r  = 0.05
nu = 1.5
D1 = 0.075
D2 = 20.0
Lx = 70.0
#
# Re-scaled parameter values needed
#
beta  = r/k1
tau   = (c+r)/k1
kappa = 0.5*(1-beta/tau)
eps   = sqrt(D1/(Lx**2*(c+r)))
D     = D2/(Lx**2*k1)
D0    = eps*D
#
# Terms x0-dependent and nonlocal coefficients
#
alpha = lambda z: exp(-nu*z)
gam   = lambda z: ((c+r)*k1**2)/(z*b**2)
#
GmsN = lambda zm,zx0,zs: cosh(sqrt(zs)*zm*zx0)*cosh(sqrt(zs)*zm*(1-zx0))
GmsD = lambda zm,zx0,zs: sqrt(zs)*zm*sinh(sqrt(zs)*zm)
Gms  = lambda zm,zx0,zs: GmsN(zm,zx0,zs)/GmsD(zm,zx0,zs)
#
Xi = lambda zga,zx0: (tau*alpha(zx0))/(36*beta**2*zga)
#
MuN = lambda zm,zga,zx0,zs: 12*Xi(zga,zx0)*Gms(zm,zx0,zs)
MuD = lambda zm,zga,zx0,zs: D0+6*Xi(zga,zx0)*Gms(zm,zx0,zs)
Mu  = lambda zm,zga,zx0,zs: MuN(zm,zga,zx0,zs)/MuD(zm,zga,zx0,zs)
#
ThetaP  = lambda lg,zx0,zm,zs,zga,eps: lg+1+zs*eps**2*zm**2-2*kappa
ThetaQ  = lambda lg,zx0,zm,zs,zga,eps: lg+1+zs*eps**2*zm**2-kappa*Mu(zm,zga,zx0,zs)
ThetaPQ = lambda lg,zx0,zm,zs,zga,eps: ThetaP(lg,zx0,zm,zs,zga,eps)/ThetaQ(lg,zx0,zm,zs,zga,eps)
Theta   = lambda lg,zx0,zm,zs,zga,eps: Mu(zm,zga,zx0,zs)*ThetaPQ(lg,zx0,zm,zs,zga,eps)
#
# Location from slow-dynamics ODE
gga = lambda z: ((0.5-z)*alpha(z))/(6*beta*nu*D0)
#
#---------------
#
# Define NLEP computation
mlow = 1.4
def NonLocalEigenvalueProblem(x0,k2,s,filename):
	# Initialise computations
#
	N   = 100
	mup = sqrt(5/4)/(sqrt(s)*eps)+20
#	
	mValues = linspace(mlow,mup,N,endpoint=True)
	mValues = mValues[::-1]
#
	numSolutions = len(mValues)
	lambdaValues = zeros(size(mValues))
#	
	fout = open(filename,'w')
	for ii in range(numSolutions):
		m = mValues[ii]
#	
		# Secondary parameters
#	
		gamma0 = gam(k2)
		Gm0    = Gms(m,x0,s)
#	
		# Lambda0 iteration
#	
		if ii == 0:
			lambda0 = 0
		else:
			lambda0 = lambdaValues[ii-1]
#		
			# Compute spatial operators
#		
			nx = 600
			x  = linspace(-Lx,Lx,nx,endpoint=True)
			hx = abs(x[1]-x[0])
#		
			# Second spatial derivative with Neumann conditions
#		
			diags        = array([-1,0,1])
			data         = ones((3,nx))
			data[0][0]   = 2
			data[1][:]   = -2
			data[-1][-1] = 2
#		
			Dxx = sparse.spdiags(data,diags,nx,nx)
			Dxx = Dxx/hx**2
#		
			# Integral weights
#		
			rho     = ones(size(x))
			rho[0]  = 1./2.
			rho[-1] = 1./2.
			rho     = rho*hx
#		
			# Compute integral of w(xi)**2
#		
			w  = 1.5/(cosh(x/2.))**2
			Iw = sum(w**2*rho)
#		
			# Compute theta
#		
			Xi0    = Xi(gamma0,x0)
			Mu0    = Mu(m,gamma0,x0,s)
			Thetah = Theta(lambda0,x0,m,s,gamma0,eps)
#										
			# Assemble eigenvalue problem
#		
			NM = Thetah/Iw*w**2*tile(w*rho,(nx,1))
			M  = Dxx+diag(2*w-1)-NM
#		
			Evals,Evec = linalg.eig(M)
			zetaR      = sort(real(Evals))
			zetaI      = sort(imag(Evals))
#		
			zetaRmax  = max(zetaR)
			izetaRmax = max(zetaI)
#		
			lambdaValues[ii] = zetaRmax-s*eps**2*m**2
#			
			print "lambda = ",lambdaValues[ii],"..."
#			
			# Save data in a readable file
#		
			fout.write(str(m))
			fout.write('\t')
			fout.write(str(lambdaValues[ii]))
			fout.write('\n')
	fout.close()
#
#---------------
#
# Compute NLEP for several s and x0, and...
#
labelbif = [0.35,0.375,0.4]
#
col = ['#339966','#3399FF','#6600CC','#FF3300']
sss = [4.5,5,5.5,6,6.5]
#
mupmax = sqrt(5/4)/(sqrt(min(sss))*eps)
#
# ... plot dispersion relation for several s, and...
#
figure(0)
for jj in range(3):
	lb      = labelbif[2]
	xnaught = lb#x0Range[lb]
	k20     = gam(gga(xnaught))#k20Range[lb]
#
	fname0 = str('nlep_s%0.3d'%jj)+'.txt'
#	
	print "Plotting dispersion relation for s = ", sss[jj], "..."
#	
	NonLocalEigenvalueProblem(xnaught,k20,sss[jj],fname0)
#	
	print "Done"
#
	nlep0 = loadtxt(fname0)
	mm0   = nlep0[:,0]
	ll0   = nlep0[:,1]
#	
	if jj == 0:
		ymax0 = max(ll0)+0.01
	else:
		ymax0 = ymax0
#
	labelname0 = str(r'$s=%0.2f$'%sss[jj])
	plot(mm0/pi,ll0,color=col[jj],lw=2,label=labelname0)
	legend(fancybox=True,shadow=True,loc=0)
#	
	# Remove .txt files (comment on if desire not)
	os.remove(fname0)
#	
title(r'$x_0 = %0.2f$'%xnaught)
mx = xlabel(r'$k$')
ly = ylabel(r'$\Re(\lambda)$')
ylim(0,ymax0)
xlim(mlow,sqrt(mupmax/pi))	
grid(True)
hold(True)
savefig('disrelS.pdf')
#
print "Done"
#
# ... plot dispersion relation for several x0
#
figure(1)
for jj in range(3):
	lb      = labelbif[jj]
	xnaught = lb
	k20     = aux(gga(xnaught))
#
	fname1 = str('nlep_x0%0.4d'%jj)+'.txt'
#	
	print "Plotting dispersion relation for x0 = ", xnaught,"..." 
#	
	NonLocalEigenvalueProblem(xnaught,k20,sss[2],fname1)
#	
	print "Done"
#
	nlep1 = loadtxt(fname1)
	mm1   = nlep1[:,0]
	ll1   = nlep1[:,1]
#	
	if jj == 0:
		ymax1 = max(ll1)+0.01
	else:
		ymax1 = ymax1
#
	labelname1 = str(r'$x_0=%0.3f$'%xnaught)
	plot(mm1/pi,ll1,color=col[jj],lw=2,label=labelname1)
	legend(fancybox=True,shadow=True,loc=0)
#	
	# Remove .txt files
	os.remove(fname1)
#	
title(r'$s = %0.2f$'%sss[2])
mx = xlabel(r'$k$')
ly = ylabel(r'$\Re(\lambda)$')
xlim(mlow-0.25,sqrt(mupmax/pi))
ylim(0,ymax1)
grid(True)
hold(True)
savefig('disrelX0.pdf')
#
#---------------
#
show()