USE rbmat
USE cbmat
USE gnuplot
! USE rtchecks
!(  k_B=1.3806503E-23 m^2 kg s^-2 K^-1
    T= 300 K
    rho=1.1174 kg/m^3
    ni=1.568E-5 m^2/s
    kTrho^-1ni^-2=1.508E-11 m
!)
kTsurhoni2d=1.508E-8
INTEGER Re
! Re=1600
ni==1/Re
ymax=50; ny=300
xmin=0.05; xmax=1; nx=50
omin=0.03; omax=0.065; nomg=20
!omin=0.03; omax=0.1; nomg=20
bmax=0.12; nbeta=10
UseContinuousAdjoint=NO

dx=(xmax-xmin)/nx
domega=(omax-omin)/nomg
!domega=omax/(nomg+0.5); omin=0.5*domega
dbeta=bmax/(nbeta+0.5); bmin=0.5*dbeta

ARRAY(-1..ny+1) OF REAL y
DO y(i)=ymax*i/ny*(0.8*i/ny+0.2) FOR ALL i

STRUCTURE[ARRAY(-2..2) OF REAL d0,d1,d2,d3,d4] derivatives(1..ny-1)
MODULE setup_derivatives
REAL M(0..4,0..4),t(0..4)
LOOP FOR iy=1 TO ny-1 WITH derivatives(iy)
  DO M(i,j)=(y(iy-2+j)-y(iy))**(4-i) FOR ALL i AND j=0 TO 4; LUdecomp M
  t=0; t(0)=24
  d4(-2+(*))=M\t
  DO M(i,j)=(5-i)*(6-i)*(7-i)*(8-i)*(y(iy-2+j)-y(iy))**(4-i) FOR ALL i AND j=0 TO 4; LUdecomp M
  DO t(i)=SUM {d4(j)*(y(iy+j)-y(iy))**(8-i)} FOR j=-2 TO 2 FOR ALL i
  d0(-2+(*))=M\t
  DO M(i,j)=(y(iy-2+j)-y(iy))**(4-i) FOR ALL i AND j=0 TO 4; LUdecomp M
  t=0; DO t(i)=SUM d0(j)*(4-i)*(3-i)*(2-i)*(y(iy+j)-y(iy))**(1-i) FOR j=-2 TO 2 FOR i=0 TO 1
  d3(-2+(*))=M\t
  t=0; DO t(i)=SUM d0(j)*(4-i)*(3-i)*(y(iy+j)-y(iy))**(2-i) FOR j=-2 TO 2 FOR i=0 TO 2
  d2(-2+(*))=M\t
  t=0; DO t(i)=SUM d0(j)*(4-i)*(y(iy+j)-y(iy))**(3-i) FOR j=-2 TO 2 FOR i=0 TO 3
  d1(-2+(*))=M\t
REPEAT
END setup_derivatives

REAL x
ARRAY(0..ny) OF REAL U0,U1,U2,V
SUBROUTINE BaseFlow()
  REAL F(0..ny)
  F(0)=0; U0(0)=0; U0(1)=0.332/SQRT(x)*y(1); U0(ny)=1
  DO
    U0(1)=U0(1)*U0(ny)**(-1.5)
    LOOP FOR i=1 TO ny-1
      d1=2/[y(i+1)-y(i-1)]; d2p=d1/[y(i+1)-y(i)]; d2m=d1/[y(i)-y(i-1)]; d20=d2p+d2m
      F(i)=F(i-1)+[U0(i-1)+U0(i)]*[y(i)-y(i-1)]/2
      U0(i+1)={d20*U0(i)-[d2m-d1*F(i)/4/x]*U0(i-1)}/[d2p+d1*F(i)/4/x]
      U1(i)=d1*[U0(i+1)-U0(i-1)]/2
      U2(i)=d2p*U0(i+1)-d20*U0(i)+d2m*U0(i-1)
    REPEAT LOOP
!    WRITE "blasius:" U0(1)/y(1),U0(ny)
  FOR 2 TIMES
  U2(0)=0; U2(ny)=0; F(ny)=F(ny-1)+[U0(ny-1)+U0(ny)]*[y(ny)-y(ny-1)]/2
  DO V(i)=(y(i)*U0(i)-F(i))/2/SQRT(x)/Re FOR i=0 TO ny
END BaseFlow

REAL omega,beta; COMPLEX alpha
ARRAY(1..ny-1,-2..2) OF COMPLEX A,dA
SUBROUTINE BuildMats()
  k2=alpha*alpha+beta*beta; k=SQRT(k2)
  LOOP FOR i=1 TO ny-1 WITH derivatives(i)
IF UseContinuousAdjoint THEN
    A(i)=[y(i+1)-y(i-1)]/2*{[I*omega-I*alpha*U0(i)]*(k2*d0-d2)+2*I*alpha*U1(i)*d1+ni*(d4-2*k2*d2+k2*k2*d0)}
    dA(i)=[y(i+1)-y(i-1)]/2*{-I*U0(i)*(k2*d0-d2)+2*I*U1(i)*d1+
      [I*omega-I*alpha*U0(i)]*2*alpha*d0+ni*(-4*alpha*d2+4*alpha*k2*d0)}
ELSE
    A(i)=[y(i+1)-y(i-1)]/2*{[I*omega-I*alpha*U0(i)]*(k2*d0-d2)-I*alpha*U2(i)*d0+ni*(d4-2*k2*d2+k2*k2*d0)}
    dA(i)=[y(i+1)-y(i-1)]/2*{-I*U0(i)*(k2*d0-d2)-I*U2(i)*d0+
      [I*omega-I*alpha*U0(i)]*2*alpha*d0+ni*(-4*alpha*d2+4*alpha*k2*d0)}
END IF
  REPEAT LOOP
  A(1,0)=A(1,0)+A(1,-2)
  dA(1,0)=dA(1,0)+dA(1,-2)
IF FALSE THEN
  A(ny-1,0)=A(ny-1,0)+A(ny-1,2)
  dA(ny-1,0)=dA(ny-1,0)+dA(ny-1,2)
ELSE
  A(ny-1,0)=A(ny-1,0)+A(ny-1,2)*EXP(-k*(y(ny+1)-y(ny-1)))+A(ny-1,1)*EXP(-k*(y(ny)-y(ny-1)))
  dA(ny-1,0)=dA(ny-1,0)+dA(ny-1,2)*EXP(-k*(y(ny+1)-y(ny-1)))+dA(ny-1,1)*EXP(-k*(y(ny)-y(ny-1)))
  A(ny-2,1)=A(ny-2,1)+A(ny-2,2)*EXP(-k*(y(ny)-y(ny-1)))
  dA(ny-2,1)=dA(ny-2,1)+dA(ny-2,2)*EXP(-k*(y(ny)-y(ny-1)))
END IF
  LUdecomp A
END BuildMats

ARRAY(-1..ny+1) OF COMPLEX zst=0,ps=0
zstar==zst(1..ny-1)
psi==ps(1..ny-1)
ARRAY(1..ny-1) OF COMPLEX zeta
SUBROUTINE eigenval()
  TYPEOF(zeta) zetanew
  DO
    BuildMats
IF UseContinuousAdjoint THEN
    psi=zeta/A; zetanew=psi*dA
    zstar=dA*zstar; zstar=A\zstar
ELSE
    psi=A\zeta; zetanew=dA*psi
    zstar=zstar*dA; zstar=zstar/A
END IF
    dalpha=-(zstar*zeta)/(zstar*zetanew)
    alpha=alpha+dalpha
    normfact=1/psi(1)
    psi=psi*normfact
    zeta=zetanew*normfact
    zstar=zstar/(zstar*zeta)
    WRITE ":" ./.
  WHILE ABS(dalpha)>1.E-8
  ps(-1)=ps(1); zst(-1)=zst(1)
  WRITE "." ./.
END eigenval

!(
v_y-I*alpha*u-I*beta*w=0
(sigma-I*alpha*U0)*u+U0_y*v-I*alpha*p-ni*lapl(u)=-I*alpha*T11+T21_y-I*beta*T31
(sigma-I*alpha*U0)*v+p_y-ni*lapl(v)=-I*alpha*T12+T22_y-I*beta*T32
(sigma-I*alpha*U0)*w-I*beta*p-ni*lapl(w)=-I*alpha*T13+T23_y-I*beta*T33

va_y=-I*alpha*ua-I*beta*wa
(sigma_0-I*alpha*U0)*ua-I*alpha*ca-ni*lapl(ua)=0
(sigma_0-I*alpha*U0)*va+U0_y*ua-ca_y-ni*lapl(va)=0
(sigma_0-I*alpha*U0)*wa-I*beta*ca-ni*lapl(wa)=0

(sigma-sigma_0)(ua*u+va*v+wa*w)=
=ua*(-I*alpha*T11+T21_y-I*beta*T31)+va*(-I*alpha*T12+T22_y-I*beta*T32)+wa*(-I*alpha*T12+T22_y-I*beta*T32)=
-I*alpha*ua*T11-ua_y*T21-I*beta*ua*T31+
-I*alpha*va*T12-va_y*T22-I*beta*va*T32+
-I*alpha*wa*T13-wa_y*T23-I*beta*wa*T33

A^2=kTmu*[4*|alpha|^2*|ua|^2+4*|va_y|^2+4*beta^2*|wa|^2+
+|ua_y+I*alpha*va|^2+|beta*ua+alpha*wa|^2+|wa_y+I*beta*va|^2]

etaa=I*beta*ua-I*alpha*wa
(alpha^2+beta^2)*ua=-I*beta*etaa+I*alpha*va_y
(alpha^2+beta^2)*wa=I*alpha*etaa+I*beta*va_y
(sigma_0-I*alpha*U0)*etaa-ni*lapl(etaa)=0
(sigma_0-I*alpha*U0)*va_y-(alpha^2+beta^2)*ca-ni*lapl(va_y)=0
(sigma_0-I*alpha*U0)*[va_yy-k^2*va]-2*I*alpha*U0_y*va_y-ni*lapl(lapl(va))=-U0_y*I*beta*etaa

A^2=kTmu*[4*|alpha|^4/|k|^4*|va_y|^2+4*|va_y|^2+4*beta^4/|k|^4*|va_y|^2+
+|alpha|^2*|va_yy/k^2+va|^2+4*beta^2*alpha^2/|k|^4*|va_y|^2+beta^2*|va_yy/k^2+va|^2]=
=kTmu*[|va_yy+k^2*va|^2*(|alpha|^2+beta^2)+4*|va_y|^2*(|alpha|^4+beta^4+|k|^4+beta^2*|alpha|^2)]/|k|^4=
=kTmu*{|va_yy|^2*(|alpha|^2+beta^2)/|k|^4+|va|^2*(|alpha|^2+beta^2)+|va_y|^2*[4*|alpha|^4+4*beta^4+4*|k|^4+4*beta^2*|alpha|^2-2*Re(k^2)*(|alpha|^2+beta^2)]/|k|^4}
!)
REAL FUNCTION SpectRec(INTEGER iy)
k2=NORM(alpha)+beta*beta
k4=NORM(alpha^2+beta^2)
WITH derivatives(iy) RESULT=NORM[SUM d2(j)*zst(iy+j) FOR j=d2.LO TO d2.HI]*k2/k4+
   k2*NORM[SUM d0(j)*zst(iy+j) FOR j=d0.LO TO d0.HI]+
   NORM[SUM d1(j)*zst(iy+j) FOR j=d1.LO TO d1.HI]*2*k2*[([IMAG(alpha)]^2+k2)/k4+4]
END SpectRec

SUBROUTINE initef
  alpha=omega/0.3
  DO zeta(i)=CRAND(); zstar(i)=CRAND() FOR  ALL i
  BaseFlow
  BuildMats
IF UseContinuousAdjoint THEN
  DO zeta=zeta/A; zstar=A\zstar FOR 2 TIMES
ELSE
  DO zeta=A\zeta; zstar=zstar/A FOR 2 TIMES
END IF
  eigenval()
END initef

!(
omega=0.05; beta=0
x=xmax
initef
COMPLEX fact
REAL N=0
WRITE x,alpha
LOOP WHILE IMAG(alpha)>0
  zstarold=zstar; alphaold=alpha
  x=x-dx
  BaseFlow(); eigenval()
  WRITE x,alpha
  fact=(zstarold*zeta)/(zstar*zetaold)*EXP[-I*(alpha+alphaold)/2/ni*dx]
  N=N+[IMAG(alpha)+IMAG(alphaold)]/2*Re*dx
  zstar=fact*zstar
REPEAT LOOP
WRITE N,LOG[ABS(fact)],0.5*LOG[0.01*0.05*0.3*Re*SUM SpectRec(iy)*[y(iy+1)-y(iy-1)]/2 FOR iy=1 TO ny-1],
      0.5*LOG(kTsurhoni2d/Re/Re)
STOP
!)
!(
Re=2100
beta=0
Nfile=CREATE("Nfile.out")
LOOP FOR omega=omin TO omax+domega/2 BY domega
  x=xmin
  initef
  REAL N=0
  LOOP WHILE x<xmax-dx/2
    alphaold=alpha
    x=x+dx
    BaseFlow(); eigenval()
    IF IMAG(alpha)>0 THEN
      IF IMAG(alphaold)<0 THEN WRITE TO Nfile SQRT(x-dx)*Re,N
      N=N+[IMAG(alpha)+IMAG(alphaold)]/2*Re*dx
      WRITE TO Nfile SQRT(x)*Re,N
    END IF
  REPEAT
  WRITE TO Nfile
REPEAT

!! N=5+6.8*(Re-1200)/1000

STOP
!)

TYPEOF(zstar) zstarold; TYPEOF(zeta) zetaold
COMPLEX FUNCTION firstappr()
  RESULT=(zstarold*zeta)/(zstar*zetaold)
END firstappr

Refile=CREATE("Refile.out")
LOOP FOR Re=1600 TO 2200 BY 100
  FILE outfile
  IF Re=1900 THEN outfile=CREATE("thermalrec.out")
  REAL vsqfin(1..ny-1),vprofile(0..ny)=0
  LOOP FOR omega=omin TO omax+domega/2 BY domega
   LOOP FOR beta=bmin TO bmax+dbeta/2 BY dbeta
    WRITE BY NAME omega,beta
    FILE xprofile
    IF Re=1900 AND ABS(omega-0.044)<1E-6 AND beta=bmin THEN xprofile=CREATE('xprofile.out')
    REAL ampl=0
    x=xmax
    initef
    WRITE x,alpha
    DO WITH derivatives(iy)
      vsqfin(iy)=NORM{[SUM d1(j)*ps(iy+j) FOR j=d1.LO TO d1.HI]/alpha}
       ! dovrebbe essere *alpha/k2
    FOR iy=1 TO ny-1
    LOOP WHILE x>xmin+dx/2
      zstarold=zstar; zetaold=zeta; alphaold=alpha
      x=x-dx
      BaseFlow(); eigenval()
      WRITE x,alpha
      fact=firstappr()*EXP[-I*(alpha+alphaold)/2/ni*dx]
      zstar=fact*zstar
      DO WITH derivatives(iy)
        IF Re=1900 AND ABS(omega-0.044)<1E-6 AND beta=bmin THEN WRITE TO xprofile x, SpectRec(10)
        ampl=ampl+SpectRec(iy)*[y(iy+1)-y(iy-1)]*dx
      FOR iy=1 TO ny-1
    REPEAT LOOP
  !  PLOT ABS(zstar(*)),zstar.LO,zstar.HI
  !  WRITE ampl
  !  READ
    DO vprofile(iy)=~+ampl*vsqfin(iy)*kTsurhoni2d/(PI*PI)/(Re*Re)*domega*dbeta FOR iy=1 TO ny-1
    IF Re=1900 THEN WRITE TO outfile omega,beta,ampl
   REPEAT LOOP
   IF Re=1900 THEN WRITE TO outfile
  REPEAT LOOP
  IF Re=1900 THEN CLOSE outfile
  WRITE TO Refile Re,MAX SQRT[vprofile(iy)] FOR iy=0 TO ny,Re*Re
  IF Re=1900 THEN
    profile=CREATE('yprofile.out')
    DO WRITE TO profile y(iy),SQRT(vprofile(iy)) FOR iy=0 TO ny
  END IF
  RANGE 0..10,0..MAX SQRT(vprofile(iy)) FOR iy=0 TO ny
  STARTLINE
  DO DRAW y(iy),SQRT(vprofile(iy)) FOR iy=0 TO ny
  SHOWGRAPH
  ! READ
REPEAT