#`help/text/ITERALL`:=TEXT(
#`HELP FOR:numerical solution of the advection partial differential`,
#`         eqaution for any one of FTCS, LAX, UPWIND, LEAPFROG`,
#`         or LAX-WENDROFF techniques`,
#` `,
#`CALLING SEQUENCE:`,
#`   ITERALL();`,
#` `,
#`SYNOPSIS:`,
#`-no parameters are required as input`,
#`-all input is prompted`,
#`-having chosen an iterative scheme this must then be entered`,
#` as a function`,
#` `,
#`EXAMPLE:`,
#`To solve the advection equation numerically using the LAX-WENDROFF`,
#`scheme, select this scheme when prompted by typing`,
#`         LAXWEND`,
#`then enter the iterative scheme at the next prompt`,
#`         (UnjP1,Unj,UnjM1,a)->Unj-`,
#`                              a*((UnjP1+Unj)/2`,
#`                                -(UnjP1+Unj)*a/2`,
#`                                -(Unj+UnjM1)/2`,
#`                                +(Unj+UnjM1)*a/2);`,
#`where UnjP1=u[n,j+1]`,
#`The numerical solution for the initial square wave evolving`,
#`according to the advection equation on the specified mesh is then`,
#`calculated and shown as a 3 dimensional plot`
#):


#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

ITERALL:=proc(a,NN,JJ,MAT,edgeS,edgeL,ans,isch)
      local j,n;

#compare FTCS(1), LAX(2), UPWIND(3), LEAPFROG(4), LAX-WENDROFF(5)

#a:=readstat(`vdt/dx__`);
#NN:=readstat(`no of time steps__`);
#JJ:=readstat(`no of space steps__`);
#MAT:=array(1..NN,1..JJ);

#initialise grid
print(`initialising`);
for n from 1 to NN do
  for j from 1 to JJ do
    MAT[n,j]:=0.0;
  od;
od;

#starting wave
#edgeS:=readstat(`start of square wave__`);
#edgeL:=readstat(`end of square well__`);

#define initial wave (square)
for j from edgeS to edgeL do
  MAT[1,j]:=1.0;
od;

#what kind of iterative scheme??????????
#ans:=readstat(`iterative scheme? FTCS, LAX, UPWIND, LEAPFROG,LAXWEND`);
#isch:=readstat(`what is this scheme?`);

if ans=1 then
for n from 2 to NN do
  for j from 2 to JJ-1 do
#   MAT[n,j]:=MAT[n-1,j]-(MAT[n-1,j+1]-MAT[n-1,j-1])*a/2;
    MAT[n,j]:=isch(MAT[n-1,j],MAT[n-1,j+1],MAT[n-1,j-1],a);
  od;
od;
fi;

if ans=2 then
for n from 2 to NN do
  for j from 2 to JJ-1 do
#   MAT[n,j]:=(MAT[n-1,j+1]+MAT[n-1,j-1])/2
#           -(MAT[n-1,j+1]-MAT[n-1,j-1])*a/2;
    MAT[n,j]:=isch(MAT[n-1,j+1],MAT[n-1,j-1],a);
  od;
od;
fi;

if ans=3 then
for n from 2 to NN do
  for j from 2 to JJ-1 do
    if a>0 then
#    MAT[n,j]:=MAT[n-1,j]-(MAT[n-1,j]-MAT[n-1,j-1])*a;
    MAT[n,j]:=isch(MAT[n-1,j],MAT[n-1,j-1],a);
   else
#    MAT[n,j]:=MAT[n-1,j]-(MAT[n-1,j+1]-MAT[n-1,j])*a;
     MAT[n,j]:=isch(MAT[n-1,j],MAT[n-1,j+1],a);
     fi;
   od;
 od;
 fi;
 
 if ans=4 then
#extra starting wave
for j from edgeS to edgeL do
  MAT[2,j]:=1;
od;

for n from 3 to NN do
  for j from 2 to JJ-1 do
#    MAT[n,j]:=MAT[n-2,j]-(MAT[n-1,j+1]-MAT[n-1,j-1])*a;
    MAT[n,j]:=isch(MAT[n-2,j],MAT[n-1,j+1],MAT[n-1,j-1],a);
  od;
od;
fi;

if ans=5 then
for n from 2 to NN do
  for j from 2 to JJ-1 do
#    MAT[n,j]:=MAT[n-1,j]-a*(1/2*(MAT[n-1,j+1]+MAT[n-1,j])
#                           -a/2*(MAT[n-1,j+1]+MAT[n-1,j])
#                           -1/2*(MAT[n-1,j]+MAT[n-1,j-1])
#                           +a/2*(MAT[n-1,j]+MAT[n-1,j-1]));
    MAT[n,j]:=isch(MAT[n-1,j+1],MAT[n-1,j],MAT[n-1,j-1],a);
  od;
od;
fi;


end:

#matrixplot(MAT,title=ans,
             
#heights=HISTOGRAM,axes=BOXED,labels=[`time`,`x`,`amplitude`]);

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

CROSS1:=proc(NN,JJ,MAT,pltvec,Nplt)
           local p,i,j,k,Q;

NP:=iquo(NN,Nplt);#print(NP);

k:=0;
#jumping thro time loop
for i from 1 by NP to NN do
k:=k+1;

L:=[];
#scanning thro space loop
for j from 1 to JJ do
l:=[[j,MAT[i,j]]];
L:=[op(L),op(l)];
od;

pltvec[k]:=plot(L,title=`waves at time intervals`);

od;

print(`plot structures in vector`);

end:

#iterative schemes for numerical solution of the advection equation
#FTCS
iFTCS:=(Unj,UnjP1,UnjM1,a)->Unj-(UnjP1-UnjM1)*a/2:
#LAX
iLAX:=(UnjP1,UnjM1,a)->(UnjP1+UnjM1)/2-(UnjP1-UnjM1)*a/2:
#UPWIND
iUPWINDp:=(Unj,UnjM1,a)->Unj-(Unj-UnjM1)*a:
iUPWINDm:=(Unj,UnjP1,a)->Unj-(UnjP1-Unj)*a:
#LEAPFROG
iLEAP:=(UnM1j,UnjP1,UnjM1,a)->UnM1j-(UnjP1-UnjM1)*a:
#LAX WENDROFF
iLWEND:=(UnjP1,Unj,UnjM1,a)->Unj-
                             a*((UnjP1+Unj)/2
                               -(UnjP1+Unj)*a/2
                               -(Unj+UnjM1)/2
                               +(Unj+UnjM1)*a/2):

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

with(plots):with(student):