USE cbmat
USE symbolic
USE gnuplot
OPENGRAPH
USE rtchecks
nx=100; ny=790
alphamin=1E-7
alphamax=10
REAL Re=400
confined=NO
external=NO
REAL y(-1..ny+1),U(-1..ny+1),nu(0..ny+1)
IF external THEN
  bf=OPEN(COMMANDLINE(1))
  DO READ FROM bf y(iy),U(iy) FOR iy=0 TO ny
  y(-1)=2*y(0)-y(1)
  U(-1)=3*U(0)-3*U(1)+U(2)
  y(ny+1)=2*y(ny)-y(ny-1)
  U(ny+1)=3*U(ny)-3*U(ny-1)+U(ny-2)
END IF
COMPLEX psi(-1..ny-1),A(-1..ny-1,-2..2)
REAL alp(0..nx); COMPLEX up(0..nx)
a0=LOG(3.109)/0.392+4.48
REAL FUNCTION f(REAL z)=LOG(z+3.109)/0.392+4.48-(a0+(0.4930-0.02450*z)*z)/(1+(0.05736+0.01101*z)*z)*EXP(-0.03385*z)
REAL FUNCTION fracc(REAL z)=IF z>3 THEN f(z) ELSE z+[f(3)-3]*(z/3)^2
LOOP FOR ix=0 TO nx
  alpha=alphamin*(alphamax/alphamin)^(ix/nx); alp(ix)=alpha
!  Re=0.5/alpha ! 2*PI/alpha
  Ly=IF 40/alpha^0.5<100 THEN 40/alpha^0.5 ELSE IF confined THEN Re ELSE 1/alpha
  psi=0
  IF NOT external THEN LOOP FOR iy=LO TO HI
!     y(iy)=Re*iy/ny
    y(iy)=Ly*[EXP(LOG[Ly*10]*iy/ny)-1]/[Ly*10-1]
!     y(iy)=Ly*[0.01*iy/ny+0.99*(iy/ny)^2]
!     U(iy)=IF y(iy)<Re THEN y(iy)-y(iy)^2/2/Re ELSE Re/2
!     U(iy)=Re/2*erf[y(iy)*SQRT(PI)/Re]
!     U(iy)=IF y(iy)<Re*PI/4 THEN Re/2*SIN(y(iy)*2/Re) ELSE Re/2
    U(iy)=fracc(y(iy)); IF confined THEN U(iy)=~+y(iy)/Re-0.57*[y(iy)/Re]^7
  REPEAT
  DO nu(iy)={IF confined THEN 1-y(iy)/Re ELSE 1}*[y(iy)-y(iy-1)]/[U(iy)-U(iy-1)]
  FOR iy=0 TO ny+1
  #define D2(p,j) {[p(j+1)-p(j)]/[y(j+1)-y(j)]-[p(j)-p(j-1)]/[y(j)-y(j-1)]}*2/[y(j+1)-y(j-1)]
  #define D1(p,j) [p(j+1)-p(j-1)]/[y(j+1)-y(j-1)]
  LOOP FOR iy=1 TO ny-1
    #define omega(j) {D2(psi,j)-alpha^2*psi(j)}
    #define nuomg(psi,n) {[nu(n)+nu(n+1)]/2*D2(psi,n)+[nu(n+1)-nu(n)]*2/[y(n+1)-y(n-1)]*D1(U,n)/[U(n)+1E-8]*psi(n)}
    #define nuomega(n) nuomg(psi,n)-[nu(n)+nu(n+1)]/2*alpha^2*psi(n)
    #define nuOmega(n) nuomg(U,n)
!(
    #define nuomega(j) [nu(j)+nu(j+1)]/2*{omega(j)-D2(U,j)/[U(j)+1E-8]*psi(j)}+ \
    {nu(j+1)*[U(j+1)-U(j)]/[y(j+1)-y(j)]-nu(j)*[U(j)-U(j-1)]/[y(j)-y(j-1)]}*2/[y(j+1)-y(j-1)]/[U(j)+1E-8]*psi(j)
!)
    eq==I*alpha*[U(iy)*omega(iy)-D2(U,iy)*psi(iy)]+ \
      D2(nuomega,iy)-alpha^2*nuomega(iy)
    DO A(iy,j)=D[eq,psi(iy+j)] FOR j=-2 TO 2
!(
    #define nuOmega(j) [nu(j)+nu(j+1)]/2*D2(U,j)*1E-8/[U(j)+1E-8]+ \
    {nu(j+1)*[U(j+1)-U(j)]/[y(j+1)-y(j)]-nu(j)*[U(j)-U(j-1)]/[y(j)-y(j-1)]}*2/[y(j+1)-y(j-1)]/[U(j)+1E-8]*U(j)
!)
    psi(iy)=-D2(nuOmega,iy)
  REPEAT
  IF confined THEN
    A(ny-1,0)=~-A(ny-1,2)*[y(ny+1)-y(ny)]/[y(ny)-y(ny-1)]
  ELSE
    A(ny-1,0)=~+A(ny-1,1)*EXP(-alpha*[y(ny)-y(ny-1)])+A(ny-1,2)*EXP(-alpha*[y(ny+1)-y(ny-1)])
    A(ny-2,1)=~+A(ny-2,2)*EXP(-alpha*[y(ny)-y(ny-1)])
  END IF
  A(0)=0; A(0,0)=1; psi(0)=-U(0)
  A(-1)=0; A(-1,0)=1; A(-1,2)=-1; psi(-1)=U(1)-U(-1)
  LUdecomp A
  psi=A\~
  up(ix)=nu(0)*[D2(psi,0)+D2(U,0)]*I/alpha
!(
  COMPLEX phi(0..ny), psi1(0..ny); phi(0)=0
  DO
    alphainf=SQRT(alpha^2+[U(i)-U(i-1)]^2/[y(i)-y(i-1)]^2*4/[U(i)+U(i-1)]^2)
    phi(i)=phi(i-1)+alphainf*[y(i)-y(i-1)]
    psi1(i)=psi(i)*EXP(phi(i))*SQRT(alphainf)
  FOR i=100 TO ny
!)
!  plot y(0..ny):psi1.REAL w l,y(0..ny):psi1.IMAG w l
!  plot y(200..ny-2):[psi(n+1).REAL-psi(n-1).REAL+U(n+1)-U(n-1)]/[y(n+1)-y(n-1)] w l
!  READ
REPEAT
set logscale xy
GNUPLOT 'set format xy "10^{%L}"'
set grid
!set yr [0.1:*]
IF confined THEN
  plot alp:ABS(REAL(up(n)))*alp(n)*Re w l, alp:ABS(IMAG(up(n)))*alp(n)*Re w l
ELSE plot alp:ABS(REAL(up(n))) w l, alp:ABS(IMAG(up(n))) w l
!(
plot alp:ABS(up(n)) w l
READ
set nologscale y
set auto y
set format y
plot alp:ARG(up(n)) w l
!)
READ
DO WRITE alp(i),up(i) FOR i=LO TO HI