;$Id: peval.pro,v 1.2 2016/09/23 15:35:49 brandenb Exp $
;
;  This is a sample routine for computing eigenvalues and eigenvectors
;  for the Rayleigh inflection point instability problem.
;  Run this directly on the command line with ".r rayleighII.pro".
;  The eigenfunctions are plotted afterwards using ".r pevec.pro"
;  Report problems to Axel Brandenburg <Axel.Brandenburg@Colorado.edu>
;
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=5,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
;  mv idl.ps ~/tex/teach/ASTR_5410/9_Richardson_crit/fig/peval_N01.ps
;  mv idl.ps ~/tex/teach/ASTR_5410/9_Richardson_crit/fig/peval_N05.ps
;  mv idl.ps ~/tex/teach/ASTR_5410/9_Richardson_crit/fig/peval_N10.ps
;  mv idl.ps ~/tex/teach/ASTR_5410/9_Richardson_crit/fig/peval_N20.ps
;
!p.charsize=1.0
!x.margin=[7.8,.5]
!y.margin=[3.2,.5]
;
N2=.05
N2=.10
N2=.20
n=60 ;(number of mesh points)
length=10.  ;(size of the domain)
;
dk=.02  ;(interval between k values)
kmax=2.1  ;(maximum wavenumber)
;
dx=length/(n+1)
x=dx*(findgen(n)+1.)-length/2.
dx2=1./dx^2
ii=complex(0.,1.)
nk=1+fix(kmax/dk)
k=dk*findgen(nk)
;
;  profiles
;
;U=(x)*1.
U=tanh(x)
;U=x<1.>(-1.)
;U=x>(-1.)<1.
U1=deriv(x,U)
U2=deriv(x,U1)
;
evalk=complexarr(n,nk)
eveck=complexarr(n,nk)
;
;  set matrix points
;
D2=fltarr(n,n)
for i=0,n-1 do D2(i,i)=-2.*dx2
for i=0,n-2 do D2(i,i+1)=dx2
for i=1,n-1 do D2(i,i-1)=dx2
;
UU=fltarr(n,n) & for i=0,n-1 do UU[i,*]=U[i]
UU2=fltarr(n,n) & for i=0,n-1 do UU2[i,i]=U2[i]
N2U2=fltarr(n,n) & for i=0,n-1 do N2U2[i,i]=N2/U[i]^2
;
for ik=0,nk-1 do begin
  k2=k(ik)^2
  S1=D2-k2*identity(n)+N2U2
  S2=D2-k2*identity(n)-N2U2
  M=UU*S1-UU2
  eval=LA_EIGENPROBLEM(M,S2,eigenvectors=evec)
  ;
  ;  sort
  ;
  good=sort(imaginary(eval))
  evalk[*,ik]=eval(good)
  eveck[*,ik]=evec[*,good[n-1]]
  ;if ik eq 20 then stop
endfor
;
!p.multi=[0,2,1]
!x.title='!8k!6'
!y.title='!6Re !8c!6'
plot,k,float(evalk[n-1,*]),yr=[-1,1]*1.05e0,/nodata
for i=0,n-1 do oplot,k,float(evalk[i,*]),ps=3
i=0   & oplot,k,abs(float(evalk[i,*])),ps=1,col=122
i=n-1 & oplot,k,abs(float(evalk[i,*])),ps=1,col=55
;
!y.title='!6Im !8c!6'
plot,k,imaginary(evalk[n-1,*]),yr=[-1,1]*max(imaginary(evalk))*1.2,/nodata
for i=0,n-1 do oplot,k,imaginary(evalk[i,*]),ps=3
i=0   & oplot,k,imaginary(evalk[i,*]),ps=1,col=122
i=n-1 & oplot,k,imaginary(evalk[i,*]),ps=1,col=55
;
!p.multi=0
print,minmax(imaginary(evalk))
!x.title='!6'
!y.title='!6'
end