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
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)=Ly*[EXP(LOG[Ly*10]*iy/ny)-1]/[Ly*10-1]
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)
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
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
REPEAT
set logscale xy
GNUPLOT 'set format xy "10^{%L}"'
set grid
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
READ
DO WRITE alp(i),up(i) FOR i=LO TO HI