USE stokes1
! USE rtchecks
FILE profile=CREATE("wrecept")
FILE spectrum=CREATE("wspectr")
ARRAY(0..nt) OF COMPLEX psiwstar,uwstar
ARRAY(0..nt) OF REAL tt
DO  tt(it)=tin+(tfin-tin)∗(it)∗(0.5∗(it)+nt0)/nt/(0.5∗nt+nt0) FOR ALL it
no=200
ARRAY(0..no) OF REAL omega
TYPEOF(zstar) zstarold
IF PlotOn THEN GRAPHICSMODE

SUBROUTINE ft(COMPLEX f^(∗),g^(∗))
LOOP FOR i=0 TO HI(omega)
    g(i)={SUM EXP[I∗omega(i)∗tt(j)]∗f(j)∗(tt(j+1)-tt(j-1)) FOR j=LO(tt)+1 TO HI(tt)-1}/2
REPEAT LOOP
END ft

DO zeta(i)=y(i)∗EXP(-y(i)∗y(i)/tfin) FOR ALL i
LOOP FOR alfa=0.3 DOWN TO 0.05 BY 0.01
  t=tfin
  sigma=(0.8∗I)∗alfa∗Re
  zstar=zeta
  eigenval
  CLEARSCREEN
  WRITE "Direct eigenfunction: ph. speed=" sigma/(I∗alfa∗Re)
  Plot(psi,8,MAXABS(psi))
  zstar1=0
  nexttime=time(NULL)+5
  LOOP FOR it=nt DOWN TO 1
    zstarold=zstar
    MB2(it)
    psiwstar(it)={zstarold(1)∗B(1,-1)+zstarold(2)∗B(2,-2)-
                 zstar(1)∗[B(1,-1)+A(1,-1)∗dt]-zstar(2)∗[B(2,-2)+A(2,-2)∗dt]}/dt
    uwstar(it)=[zstarold(1)∗B(1,-2)-zstar(1)∗(B(1,-2)+dt∗A(1,-2))]∗(y(-1)-y(1))/dt
    IF it MOD(nt DIV 100) = 0 THEN
      WRITE TO profile alfa,t+dt,Energy(zstarnorm),ABS(psiwstar(it)),
               ARG(psiwstar(it)),ABS(uwstar(it)),ARG(uwstar(it))
    END IF
  REPEAT LOOP
  WRITE TO profile
  COMPLEX psiwstarft(0..no),uwstarft(0..no)
  DO omega(i)=(0.5+0.5∗i/no)∗alfa∗Re FOR ALL i
  ft(psiwstar,psiwstarft)
  ft(uwstar,uwstarft)
  DO WRITE TO spectrum alfa,omega(i),ABS(psiwstarft(i)),ARG(psiwstarft(i)),
              ABS(uwstarft(i)),ARG(uwstarft(i))
  FOR ALL i
  WRITE TO spectrum
REPEAT LOOP

IF PlotOn THEN TEXTMODE