! PROGRAM BAFIT (Back And Forth In Time)
USE cbmat
USE graphics
! USE rtchecks
vsq_wsq=0; usq=1
energyin=vsq_wsq
Stewartson=FALSE
fixedbeta=FALSE
wrc=FALSE
FILE INPUT=OPEN("bafit.in")
plot=BOOLEAN FROM INPUT; G=REAL FROM INPUT; nmodes=INTEGER FROM INPUT
xin=REAL FROM INPUT; dx=REAL FROM INPUT; xfin=REAL FROM INPUT
ny=INTEGER FROM INPUT; n0=INTEGER FROM INPUT; ymax=REAL FROM INPUT
IF plot THEN GRAPHICSMODE
REAL y(-1..ny)
DO y(i)=ymax∗i∗(0.5∗i+n0)/ny/(0.5∗ny+n0) FOR i=-1 TO ny
wpos=0; ppos=1; upos=2; vpos=3;
Uw=wpos; UvpVu=ppos; Uu=upos; u=vpos
REAL omega=0, beta=0

REAL FUNCTION Penergy(COMPLEX f^(∗)) = SUM NORM[f(4∗i+u)]∗[y(i)-y(i-1)] FOR i=1 TO ny-1

REAL FUNCTION Tenergy(COMPLEX f^(∗))
  RESULT=0 ! [f(4∗(ny-2)+UvpVu).REAL∗f(4∗(ny-2)+UvpVu).REAL+f(4∗(ny-2)+UvpVu).IMAG∗f(4∗(ny-2)+UvpVu).IMAG]/beta
  DO RESULT=RESULT+NORM[f(4∗i+Uw)]∗[y(i)-y(i-1)]+
            NORM[f(4∗i+UvpVu)]∗[y(i+1)-y(i-1)]/2 FOR i=1 TO ny-2
END Tenergy

ARRAY(0..4∗ny-1,-7..7) OF COMPLEX A,B,T
ARRAY(0..4∗ny-1) OF COMPLEX f,fold,fstar,fstar1
ARRAY(0..ny) OF REAL F,U,V

COMPLEX FUNCTION sclrprod = SUM fstar(i)∗f(i)+fstar1(i)∗(f(i)-fold(i)) FOR ALL i

SUBROUTINE BaseFlow(REAL x)
  F(0)=0; U(0)=0; U(1)=0.332/SQRT(x)∗y(1); U(ny)=1
  DO
    U(1)=U(1)∗U(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)+[U(i-1)+U(i)]∗[y(i)-y(i-1)]/2
      U(i+1)={d20∗U(i)-[d2m-d1∗F(i)/4/x]∗U(i-1)}/[d2p+d1∗F(i)/4/x]
    REPEAT LOOP
    !    WRITE "blasius:" U(1)/y(1),U(ny)
  FOR 2 TIMES
  ! V=(y∗U-F)/2./x
  F(ny)=F(ny-1)+[U(ny-1)+U(ny)]∗[y(ny)-y(ny-1)]/2
  DO V(i)=MIN{[y(i)∗U(i)-F(i)]/2/x,0.8605/SQRT(x)} FOR i=0 TO ny
  !  NEWGRAPH; MOVE(0,0); DO DRAW(y(i),U(i)) FOR i=1 TO ny-1; SHOWGRAPH
END BaseFlow

SUBROUTINE BuildMats
  b2=beta∗beta
  A=0; B=0
  A(0+wpos,0)=(n0+0.25)/2/n0
  A(0+wpos,4)=(n0-0.25)/2/n0
  A(0+ppos,0)=1; A(0+ppos,4)=-1
  A(0+upos,0)=1
  A(0+vpos,0)=1
  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
    !    d1m=1/[y(i)-y(i-1)]; d2mp=2∗d1/[y(i+1)-y(i-1)]; d2mm=2∗d1/[y(i)-y(i-2)]; d2m0=d2mp+d2mm
    !!!!!!! correzione del 18-10-97    
    d1m=1/[y(i)-y(i-1)]; d2mp=d1m∗d1; d2mm=2∗d1m/[y(i)-y(i-2)]; d2m0=d2mp+d2mm
    eq=4∗i
    !  i∗omega∗w+U∗wx+V∗wy-i∗beta∗p-wyy+b2∗w=0
    !  i∗omega∗w+(U∗w)x+(V∗w)y-i∗beta∗p-wyy+b2∗w=0
    A(eq+wpos,-4)=-d1m∗V(i-1)/2-d2mm
    A(eq+wpos,0)=d1m∗(V(i)-V(i-1))/2+d2m0+b2+I∗omega
    A(eq+wpos,4)=d1m∗V(i)/2-d2mp
    A(eq+wpos,ppos-wpos)=-I∗beta
    B(eq+wpos,0)=(U(i)+U(i-1))/2
    !  i∗omega∗v+U∗vx+V∗vy+Vx∗u+Vy∗v+2∗G∗G∗U∗u+py-vyy+b2∗v=0
    !  i∗omega∗v+(U∗v+V∗u)x+(2∗V∗v)y-i∗beta∗V∗w+2∗G∗G∗U∗u+py-vyy+b2∗v=0
    A(eq+ppos,vpos-ppos-4)=-d1∗V(i-1)-d2m
    A(eq+ppos,vpos-ppos)=d20+b2+I∗omega
    A(eq+ppos,vpos-ppos+4)=d1∗V(i+1)-d2p
    A(eq+ppos,upos-ppos)=2∗G∗G∗U(i)
    A(eq+ppos,wpos-ppos)=-I∗beta∗V(i)/2
    A(eq+ppos,wpos+4-ppos)=-I∗beta∗V(i)/2
    A(eq+ppos,0)=-d1
    A(eq+ppos,4)=d1
    B(eq+ppos,vpos-ppos)=U(i)
    B(eq+ppos,upos-ppos)=V(i)
    !  i∗omega∗u+(U∗u)x+V∗uy+Uy∗v-uyy+b2∗u=0
    A(eq+upos,-4)=-d1∗V(i)/2-d2m
    A(eq+upos,0)=d20+b2+I∗omega
    A(eq+upos,4)=d1∗V(i)/2-d2p
    A(eq+upos,vpos-upos)=d1∗(U(i+1)-U(i-1))/2
    B(eq+upos,0)=U(i)
    !  ux+vy-i∗beta∗w=0
    A(eq+vpos,wpos-vpos)=-I∗beta
    A(eq+vpos,-4)=-d1m
    A(eq+vpos,0)=d1m
    B(eq+vpos,upos-vpos)=0.5
    B(eq+vpos,upos-vpos-4)=0.5
  REPEAT LOOP
END BuildMats

SUBROUTINE WKB
  TYPEOF(f) f2old,f1,v1
  COMPLEX sigma=0
  REAL x=xin
  LOOP WHILE x<xfin-0.000001
    x=x+dx
    f2old=fold; fold=f
    BaseFlow(x)
    BuildMats
    DO
      T=A+sigma∗B
      LUdecomp T
      f1=T\f; f1=B∗f1
      v1=fstar∗B; v1=v1/T
      dsigma=-(v1∗f)/(v1∗f1)
      sigma=sigma+dsigma
      f=f1∗(f1∗fold)/ABS(f1∗fold)/SQRT[Penergy(f1)]
      fstar=v1/(v1∗f)
      WRITE ":" ./.
    WHILE ABS(dsigma)>1.E-5
    c1=SUM fstar(i)∗[1.5∗f(i)-2∗fold(i)+0.5∗f2old(i)]/dx FOR ALL i
    WRITE x,x∗REAL(sigma),x∗REAL(sigma-c1)
  REPEAT LOOP
END WKB

SUBROUTINE Plot(INTEGER v; REAL h,stagger,height)
  RANGE(0,8,0,height)
  DO LINE(i,0,i,height) FOR i=0 TO 8
  DO LINE(0,yy,8,yy) FOR yy=0 TO height+0.001 BY 0.2∗height
  MOVE 0,0
  DO DRAW([y(i)+stagger∗(y(i-1)-y(i))]/h,ABS(f(4∗i+v))) FOR i=1 TO ny-1
  SHOWGRAPH
END Plot

SUBROUTINE MarchForward
  ARRAY(0..4∗ny-1) OF COMPLEX f2old
  REAL x=xin
  LOOP WHILE x<xfin-0.000001
    x=x+dx
    f2old=fold; fold=f
    BaseFlow(x)
    BuildMats
    T=A+1.5/dx∗B
    DO f(i)=(2∗fold(i)-0.5∗f2old(i))/dx FOR ALL i
    LUdecomp T; f=T\f;
    f=B∗f
    !  NEWGRAPH; Plot(u,sqrt(x),0.5,sqrt(Penergy(f)/sqrt(x)))
    !  Plot(Uw,sqrt(x),0.5,sqrt(Penergy(f)/sqrt(x)))
    !  sigma={1.5+[0.5∗Penergy(f2old)-2∗Penergy(fold)]/Penergy(f)}/2/dx
    !  WRITE x,Penergy(f),x∗REAL(sigma)
  REPEAT LOOP
END MarchForward

SUBROUTINE MarchBackward
  REAL x=xfin
  LOOP WHILE x>xin+0.000001
    BaseFlow(x)
    BuildMats
    T=A+1.5/dx∗B
    fstar1=fstar+fstar1
    fstar=fstar1-fstar
    fstar1=fstar1∗B
    LUdecomp T; fstar1=fstar1/T
    fstar1=fstar1∗1.5/dx
    fstar=fstar1-fstar
    fstar1=fstar1/3
    x=x-dx
    IF wrc THEN
      DO f(4∗i+u)=CONJG[fstar(4∗i+u)+(fstar(4∗i+Uu)+fstar(4∗i+Uu+4))/2]/[y(i)-y(i-1)]
        f(4∗i+Uu)=CONJG[(fstar(4∗i+u)+fstar(4∗i+u-4))/2+fstar(4∗i+Uu)]∗2./[y(i+1)-y(i-1)]
      FOR i=1 TO ny-1
      WRITE x,Penergy(f)
    END IF  
    !  WRITE x,Tenergy(fstar)
  REPEAT LOOP
END MarchBackward

COMPLEX Mode(1..nmodes,0..ny-1)

REAL FUNCTION MaxGain(INTEGER mode)
  TYPEOF(f) fprev
  BaseFlow(xfin)
  BuildMats
  f=0; DO f(4∗i+u)=[y(i)+y(i-1)]∗[U(i)-U(i-1)]/[y(i)-y(i-1)] FOR i=1 TO ny-1
  RESULT=0
  DO
    LOOP FOR nmode=1 TO mode-1 
      COMPLEX pr=0; DO pr=pr+Mode(nmode,i)|f(4∗i+u)∗[y(i)-y(i-1)]  FOR i=1 TO ny-1
      DO f(4∗i+u)=f(4∗i+u)-pr∗Mode(nmode,i) FOR i=1 TO ny-1
    REPEAT LOOP
    prevRESULT=RESULT
    f=f/sqrt(Penergy(f))
    IF plot THEN  
      NEWGRAPH
      WRITE "final u (beta=" beta, "omega="omega, "gain: " RESULT ")"
      Plot(u,1,0.5,1)
      up==f(u+4∗(∗))
      MOVE 0,0; DO DRAW y(i),ABS{[(up(i+1)-up(i))/(y(i+1)-y(i))-(up(i)-up(i-1))/(y(i)-y(i-1))]∗2/(y(i+1)-y(i))} FOR i=2 TO ny-2
      SHOWGRAPH
    END IF  
    fprev=f
    fstar=0; fstar1=0; DO fstar(4∗i+u)=CONJG[f(4∗i+u)]∗[y(i)-y(i-1)] FOR i=1 TO ny-1
    MarchBackward
    f=0
    IF energyin=vsq_wsq THEN  
      DO f(4∗i+UvpVu)=CONJG[fstar(4∗i+UvpVu)]∗2/[y(i+1)-y(i-1)]
        f(4∗i+Uw)=-I∗[f(4∗i+UvpVu)-f(4∗i+UvpVu-4)]/[y(i)-y(i-1)]/beta 
      FOR i=1 TO ny-1
      f=f/sqrt(Tenergy(f))
      IF plot THEN  
        NEWGRAPH
        WRITE "initial v,w (beta=" beta ")"
        Plot(UvpVu,1,0,1)
        Plot(Uw,1,0.5,1)
      END IF
    ELSE IF energyin=usq THEN  
      DO f(4∗i+u)=CONJG[fstar(4∗i+u)+(fstar(4∗i+Uu)+fstar(4∗i+Uu+4))/2]/[y(i)-y(i-1)]
        f(4∗i+Uu)=CONJG[(fstar(4∗i+u)+fstar(4∗i+u-4))/2+fstar(4∗i+Uu)]∗2./[y(i+1)-y(i-1)]
      FOR i=1 TO ny-1
      f=f/sqrt(Penergy(f))
      IF plot THEN  
        NEWGRAPH
        WRITE "initial u"
        Plot(Uu,1,0.5,1)
      END IF
    END IF
    DO fold(i)=f(i) FOR ALL i
    MarchForward
    RESULT=Penergy(f)
    !  RESULT=SUM REAL{f(4∗i+u)|fprev(4∗i+u)}∗[y(i)-y(i-1)] FOR i=1 TO ny-1
    !  WRITE RESULT ./.
  UNTIL ABS(1-prevRESULT/RESULT)<1.E-5
  f=f/sqrt(Penergy(f))
  DO Mode(mode,i)=f(4∗i+u) FOR i=0 TO ny-1
END MaxGain

IF wrc THEN
  TYPEOF(f) f1,v1
  COMPLEX sigma
  beta=REAL FROM PROMPT("beta?")
  omega=REAL FROM PROMPT("omega?")
  BaseFlow(xfin)
  BuildMats
  f=0; DO f(4∗i+u)=[y(i)+y(i-1)]∗[U(i)-U(i-1)]/[y(i)-y(i-1)] FOR i=1 TO ny-1
  fstar=0; fstar1=0; DO fstar(4∗i+u)=CONJG[f(4∗i+u)]∗[y(i)-y(i-1)] FOR i=1 TO ny-1
  sigma=I∗omega/0.3
  DO
    T=A+sigma∗B
    LUdecomp T
    f1=T\f; f1=B∗f1
    v1=fstar∗B; v1=v1/T
    dsigma=-(v1∗f)/(v1∗f1)
    sigma=sigma+dsigma
    f=f1/SQRT[Penergy(f1)]
    fstar=v1/(v1∗f)
    WRITE TO stderr ":" ./.
  WHILE ABS(dsigma)>1.E-5
  MarchBackward
  DO f(4∗i+u)=CONJG[fstar(4∗i+u)+(fstar(4∗i+Uu)+fstar(4∗i+Uu+4))/2]/[y(i)-y(i-1)]
    f(4∗i+Uu)=CONJG[(fstar(4∗i+u)+fstar(4∗i+u-4))/2+fstar(4∗i+Uu)]∗2./[y(i+1)-y(i-1)]
  FOR i=1 TO ny-1
  IF plot THEN  
    NEWGRAPH
    WRITE "initial u: gain", 1/Penergy(f)
    f=f/sqrt(Penergy(f))
    Plot(Uu,1,0.5,1)
  END IF
ELSE
  
  FILE Spectrum=CREATE("spectrum")
  FILE Profile=CREATE("profile")
  IF Stewartson THEN
    BaseFlow(1)
    f=0; DO f(4∗i+u)=[y(i)+y(i-1)]∗[U(i)-U(i-1)]/[y(i)-y(i-1)] FOR i=1 TO ny-1
    f=f/sqrt(Penergy(f))
    DO
      WRITE TO Profile y(i),ABS([f(4∗i+u)+f(4∗i+4+u)]/2),
      ABS[fstar(4∗i+UvpVu)],ABS([fstar(4∗i+Uw)+fstar(4∗i+4+Uw)]/2)
    FOR i=1 TO ny-2
    STOP
  END IF
  omega=-1
  beta=0.45
  DO
    omega=omega+1
    DO
      beta=beta+0.1
      IF fixedbeta THEN
        beta=REAL FROM PROMPT("beta?")
        omega=REAL FROM PROMPT("omega?")
      END IF
      g1=MaxGain(1)
      f=f/sqrt(Penergy(f))
      DO fstar(4∗i+UvpVu)=CONJG[fstar(4∗i+UvpVu)]∗2/[y(i+1)-y(i-1)]
        fstar(4∗i+Uw)=-I∗[fstar(4∗i+UvpVu)-fstar(4∗i+UvpVu-4)]/[y(i)-y(i-1)]/beta
      FOR i=1 TO ny-1
      fstar=fstar/sqrt(Tenergy(fstar))
      !(
        WRITE TO Profile 0,0,0,0
        DO
          WRITE TO Profile y(i),ABS([f(4∗i+u)+f(4∗i+4+u)]/2),
          ABS[fstar(4∗i+UvpVu)],ABS([fstar(4∗i+Uw)+fstar(4∗i+4+Uw)]/2)
        FOR i=1 TO ny-2
        DO WRITE TO Profile y(i),0,0,0 FOR i=ny-1 TO ny
        WRITE TO Profile
        !)  
      WRITE TO Spectrum omega:15.8, beta:15.8, g1:15.8 ./.
      LOOP FOR i=2 TO nmodes; WRITE TO Spectrum ,MaxGain(i):15.8 ./.; REPEAT LOOP
      WRITE TO Spectrum
      IF fixedbeta THEN STOP
    WHILE beta<1-1.E-6
  WHILE omega<10-1.E-6
  
END IF