#include <time.h>
USE cbmat
USE graphics
FILE INPUT=OPEN("stokes1.in")
PlotOn=BOOLEAN FROM INPUT
provacostante=FALSE
Re=REAL FROM INPUT
REAL tin=REAL FROM INPUT; nt=INTEGER FROM INPUT; nt0=INTEGER FROM INPUT; REAL tfin=REAL FROM INPUT
ny=INTEGER FROM INPUT; ny0=INTEGER FROM INPUT; ymax=REAL FROM INPUT
long nexttime
REAL t,dt
REAL y(-1..ny+1)
DO y(i)=ymax*i*(0.5*i+ny0)/ny/(0.5*ny+ny0) FOR i=-1 TO ny+1
REAL alfa
COMPLEX sigma
REAL FUNCTION Energy(COMPLEX f^(*))
a2=alfa*alfa
RESULT=NORM(f(1))/y(1)+NORM(f(ny-1))/[y(ny)-y(ny-1)]
DO RESULT=RESULT+NORM(f(i)-f(i-1))/[y(i)-y(i-1)] FOR i=2 TO ny-1
DO RESULT=RESULT+a2*NORM(f(i))*[y(i+1)-y(i-1)]/2 FOR i=1 TO ny-1
END Energy
ARRAY(1..ny-1,-2..2) OF COMPLEX A,B,T
ARRAY(1..ny-1) OF COMPLEX psi,zeta,zetaold,zstar,zstar1,zstarnorm
ARRAY(0..ny) OF REAL U,U2
COMPLEX FUNCTION sclrprod()=SUM zstar(i)*zeta(i)+zstar1(i)*(zeta(i)-zetaold(i)) FOR ALL i
pi2=3.544907702
SUBROUTINE BaseFlow()
DO U(i)=erfc(y(i)/2/SQRT(t)); U2(i)=y(i)*EXP[-y(i)*y(i)/4/t]/pi2/t**1.5 FOR ALL i
END BaseFlow
SUBROUTINE BuildMats()
a2=alfa*alfa
A=0; B=0
LOOP FOR i=1 TO ny-1
d1p=2/(y(i+2)-y(i)); d2pp1=d1p/(y(i+2)-y(i+1)); d2pm1=d1p/(y(i+1)-y(i)); d2p0=-d2pm1-d2pp1
d1m=2/(y(i)-y(i-2)); d2mp1=d1m/(y(i)-y(i-1)); d2mm1=d1m/(y(i-1)-y(i-2)); d2m0=-d2mm1-d2mp1
d1=2/(y(i+1)-y(i-1))
REAL d2p1=d1/(y(i+1)-y(i)), d2m1=d1/(y(i)-y(i-1)), d20=-d2m1-d2p1
d4p2=d2p1*d2pp1
d4p1=d2p1*d2p0+d20*d2p1
d40=d2p1*d2pm1+d20*d20+d2m1*d2mp1
d4m1=d2m1*d2m0+d20*d2m1
d4m2=d2m1*d2mm1
d00=2/3; d0p1=-d2p1/3/d20; d0m1=-d2m1/3/d20
d2p2=-d4p2/6/d20
d2p1=d2p1-d4p1/6/d20
d2m2=-d4m2/6/d20
d2m1=d2m1-d4m1/6/d20
d20=d20-d40/6/d20
B(i,2)=-d2p2
B(i,1)=-d2p1+a2*d0p1
B(i,0)=-d20+a2*d00
B(i,-1)=-d2m1+a2*d0m1
B(i,-2)=-d2m2
A(i,2)=I*alfa*Re*U(i)*d2p2+d4p2-2*a2*d2p2
A(i,1)=I*alfa*Re*[U(i)*(d2p1-a2*d0p1)-U2(i)*d0p1]+d4p1-2*a2*d2p1+a2*a2*d0p1
A(i,0)=I*alfa*Re*[U(i)*(d20-a2*d00)-U2(i)*d00]+d40-2*a2*d20+a2*a2*d00
A(i,-1)=I*alfa*Re*[U(i)*(d2m1-a2*d0m1)-U2(i)*d0m1]+d4m1-2*a2*d2m1+a2*a2*d0m1
A(i,-2)=I*alfa*Re*U(i)*d2m2+d4m2-2*a2*d2m2
REPEAT LOOP
A(1,0)=A(1,0)+A(1,-2)
B(1,0)=B(1,0)+B(1,-2)
A(ny-1,0)=A(ny-1,0)+A(ny-1,2)*EXP(-alfa*(y(ny+1)-y(ny-1)))+A(ny-1,1)*EXP(-alfa*(y(ny)-y(ny-1)))
B(ny-1,0)=B(ny-1,0)+B(ny-1,2)*EXP(-alfa*(y(ny+1)-y(ny-1)))+B(ny-1,1)*EXP(-alfa*(y(ny)-y(ny-1)))
A(ny-2,1)=A(ny-2,1)+A(ny-2,2)*EXP(-alfa*(y(ny)-y(ny-1)))
B(ny-2,1)=B(ny-2,1)+B(ny-2,2)*EXP(-alfa*(y(ny)-y(ny-1)))
END BuildMats
SUBROUTINE Plot(COMPLEX f^(*); REAL h,height)
RANGE(0,h,0,height)
DO LINE(i,0,i,height) FOR i=0 TO FLOOR(h)
DO LINE(0,yy,h,yy) FOR yy=0 TO height+0.001 BY 0.2*height
STARTLINE; DO DRAW(y(i),ABS(f(i))) FOR i=0 TO ny-1
STARTLINE; DO DRAW(y(i),height*U(i)) FOR i=0 TO ny-1
WRITE
END Plot
SUBROUTINE MarchForward1
t=tin; INTEGER it=0
LOOP WHILE it<nt
it=it+1
dt=tin+(tfin-tin)*it*(0.5*it+nt0)/nt/(0.5*nt+nt0)-t
t=t+dt
BaseFlow(); BuildMats()
T=B+dt*A; LUdecomp T
psi=T\zeta; zeta=B*psi
REPEAT LOOP
END MarchForward1
SUBROUTINE MarchForward2
TYPEOF(zeta) zeta2old
t=tin; INTEGER it=0
nexttime=time(NULL)+2
LOOP WHILE it<nt
it=it+1
dt=tin+(tfin-tin)*it*(0.5*it+nt0)/nt/(0.5*nt+nt0)-t
t=t+dt
zeta2old=zetaold; zetaold=zeta
BaseFlow(); BuildMats()
T=A+(1.5/dt)*B; LUdecomp T
DO zeta(i)=(2*zetaold(i)-0.5*zeta2old(i))/dt FOR i=zeta.LO TO zeta.HI
psi=T\zeta; zeta=B*psi
IF PlotOn THEN
IF time(NULL)>=nexttime THEN
REAL psimax=ABS(psi(1))
DO psimax=MAX(psimax,ABS(psi(i))) FOR i=2 TO ny-1
CLEARSCREEN
sigma=0
DO sigma=sigma+CONJG[psi(i)]*[1.5*zeta(i)-2*zetaold(i)+0.5*zeta2old(i)]*[y(i+1)-y(i-1)]/2 FOR ALL i
sigma=sigma/Energy(psi)/dt
WRITE "time="t, "Energy="Energy(psi), "sigma="sigma
WRITE "ph.speed=" sigma/(I*alfa*Re)
Plot(psi,8,psimax)
nexttime=time(NULL)+2
END IF
END IF
REPEAT LOOP
END MarchForward2
SUBROUTINE MarchBackward1
t=tfin; INTEGER it=nt
LOOP WHILE it>0
dt=t-tin-(tfin-tin)*(it-1)*(0.5*(it-1)+nt0)/nt/(0.5*nt+nt0)
BaseFlow(); BuildMats()
T=B+dt*A; LUdecomp T
zstar=zstar*B; zstar=zstar/T
it=it-1
t=t-dt
REPEAT LOOP
END MarchBackward1
SUBROUTINE MB2(INTEGER it)
t=tin+(tfin-tin)*(it)*(0.5*(it)+nt0)/nt/(0.5*nt+nt0)
dt=t-tin-(tfin-tin)*(it-1)*(0.5*(it-1)+nt0)/nt/(0.5*nt+nt0)
BaseFlow(); BuildMats()
T=A+(1.5/dt)*B; LUdecomp T
zstar1=zstar+zstar1
zstar=zstar1-zstar
zstar1=zstar1*B; zstar1=zstar1/T
zstar1=zstar1*1.5/dt
zstar=zstar1-zstar
zstar1=zstar1/3
DO zstarnorm(i)=zstar(i)/(y(i+1)-y(i-1)) FOR i=1 TO ny-1
IF PlotOn THEN
IF time(NULL)>=nexttime THEN
REAL zstarmax=ABS(zstarnorm(1))
DO zstarmax=MAX(zstarmax,ABS(zstarnorm(i))) FOR i=2 TO ny-1
CLEARSCREEN
WRITE "alfa="alfa, "t="t, "direction: backward", "Norm="Energy(zstarnorm)
Plot(zstarnorm,8,zstarmax)
nexttime=time(NULL)+2
END IF
END IF
END MB2
SUBROUTINE MarchBackward2
nexttime=time(NULL)+2
LOOP FOR it=nt DOWN TO 1
MB2(it)
REPEAT LOOP
END MarchBackward2
SUBROUTINE eigenval()
TYPEOF(zeta) zeta1
BaseFlow; BuildMats
DO
T=A+sigma*B; LUdecomp T
psi=T\zeta; zeta1=B*psi
zstar=zstar*B; zstar=zstar/T
dsigma=-(zstar*zeta)/(zstar*zeta1)
sigma=sigma+dsigma
normfact=ABS[zeta1(1)]/zeta1(1)/SQRT[Energy(psi)]
zeta=zeta1*normfact
psi=psi*normfact
zstar=zstar/(zstar*zeta)
WRITE ":" ./.
WHILE ABS(dsigma)>1.E-8
WRITE "." ./.
END eigenval