x_message(['Integrazione problema n.2  con';...
           'schemi alle differenze finite.';' ';...
           ' '; ...
           'Ver 2.0  (27/1/2004)'])


function y=ge(x,t)
 for ii=1:N,
 y(ii)=0.;
 if (x(ii)>lambda*t & x(ii)<lambda*t+0.1) then 
    y(ii)=(sin(10*pi*(x(ii)-lambda*t))),
 end,
 end,
endfunction;

cont0="Y";
while cont0=="Y",


pi=atan (1)*4;

txt=['lambda'];
sig=x_mdialog('Parametri equazione di convezione',txt,['1.0']);
lambda=evstr(sig(1));
//mu=evstr(sig(2));
//lambda=1;
mu=0;


N=101;
dt=0.01;
itermax=70;
Resmin=1e-5;


cambia=2;
while cambia==2, 
sN=msprintf('%i', N);
sdt=msprintf('%f', dt);
sitermax=msprintf('%i', itermax);
sResmin=msprintf('%e', Resmin);
txt=['N. nodi';'dt    ';'N. max iter';'Res. min'];
//sig=x_mdialog('Parametri  integrazione ',txt,['101';'0.01';'70';'1e-5']);
sig=x_mdialog('Parametri  integrazione ',txt,[sN; sdt; sitermax; sResmin]);
N=evstr(sig(1));
dt=evstr(sig(2));
itermax=evstr(sig(3));
Resmin=evstr(sig(4));

//N=101;
dx=1/(N-1);
c=lambda*dt/dx;
s=mu*dt/dx^2;


tschema=['FTCS (esplicito,  I ord. tempo, II ord. spazio)';...
         'BTCS (implicito,  I ord. tempo, II ord. spazio)';...
         'CN   (implicito, II ord. tempo, II ord. spazio)'];
nschema=x_choose(tschema,['Scegli lo schema alle differenze finite']);


testo(1) =strcat( ['Riassunto dei parametri equazione ed  integrazione:']);
testo(2)=strcat( ['------------------------------------------------------------']);
testo(3) =strcat( ['Parametri equazione ----->  lambda:' string(msprintf('%2.5f',lambda   ))']);
testo(4) =strcat( ['                                mu:' string(msprintf('%2.5f',mu       ))']);
testo(5) =strcat( ['Parametri integrazione --> N. nodi:' string(msprintf('%2.5f',N   ))']);
testo(6) =strcat( ['                                dx:' string(msprintf('%2.5f',dx  ))']);
testo(7) =strcat( ['                                dt:' string(msprintf('%2.5f',dt  ))']);
testo(8) =strcat( ['                                 c:' string(msprintf('%2.5f',c   ))']);
testo(9) =strcat( ['                                 s:' string(msprintf('%2.5f',s   ))']);
testo(10) =strcat( ['                       N. iter max:' string(msprintf('%2.5f',itermax))']);
testo(11)=strcat( ['                         Res. min.:' string(msprintf('%2.5f',Resmin))']);
testo(12)=strcat( ['            Schema di integrazione:' string(msprintf('%s',tschema(nschema) ))']);
testo(13)=strcat( ['------------------------------------------------------------']);
testo(14)=strcat( ['Vuoi cambiare i parametri?']);
cambia= x_message(testo,['No' 'Yes']);
end,


X=[0 : dx : 1.]';
t=0.;
gn=ge(X,t);
gnp1=zeros(1,N)';


iter=2000;


if nschema==1 then
  //FTCS
  A=zeros(N,N);
  A(1,1)=1.;
  for i=2:N-1
    A(i,i+1)=s-c/2;
    A(i,i)=1-2*s;
    A(i,i-1)=(c/2+s);
  end
  A(N,N)=1.;
  L=A;
elseif nschema==2 then
  //BTCS
  A=zeros(N,N);
  A(1,1)=1.;
  for i=2:N-1
    A(i,i+1)=c/2-s;
    A(i,i)=1+2*s;
    A(i,i-1)=-(c/2+s);
  end
  A(N,N)=1.;
  L=inv(A);
else
  //CN
  A=zeros(N,N);
  B=zeros(N,N);
  A(1,1)=1.;
  B(1,1)=1.;
  for i=2:N-1
    A(i,i+1)=1/2*(c/2-s);
    B(i,i+1)=1/2*(s-c/2);
    A(i,i)=1+s;
    B(i,i)=1-s;
    A(i,i-1)=-1/2*(c/2+s);
    B(i,i-1)=1/2*(c/2+s);
  end
  A(N,N)=1.;
  B(N,N)=1.;
  L=inv(A)*B;
end,




xset("window",0)
xtitle('Evoluzione della soluzione', 'x', 'g');
//xset("window",1)
xset("window",0)
xselect();
//kp=xget("pixmap");xset("pixmap",1);
//xset("pixmap",kp);
//kp=xget("pixmap");xset("pixmap",1);

xselect();
xset("pixmap",1);
plot2d(X,gn,-1,rect=[0.,0.,1.,1.]);

xset("mark",-1,2);
t=0.;
res=1e+39;
gnp1max=0.;
nn=1;
while nn <= itermax & res > Resmin & gnp1max <= 1.2 , 
  gnp1=L*gn;
  gnp1(1)=0.;
  gnp1(N)=0.;
  gnp1max=max(gnp1);
  if gnp1max > 1.2 then
   x_message(['Attenzione   g>1' ; 'Schema di integrazione instabile ']);
   break;
  end;
//  xset("window",0)
  xset("wwpc");
  xtape("clear",0); 
  plot2d(X,gnp1,-1,rect=[0.,0.,1.,1.]);
//  xset("window",1);
//  res=sum(abs((gnp1-gn))/dt)/N;
//  plot2d(nn,res,-1,logflag="nl", rect=[0,Resmin,itermax,10.]);
  t=t+dt;
//  ge=sin(pi*X)*exp(-pi^2*t);
  plot2d(X,ge(X,t),2,rect=[0.,0.,1.,1.]);
  xset("wshow");
  gn=gnp1;
  nn=nn+1;
end
xselect();


//xset("window",0);
//ge=sin(pi*X)*exp(-pi^2*t);
//plot2d(X,ge,2,rect=[0.,0.,1.,1.]);
//xset("wshow");

//errore numerico
//xset("window",2);
//plot2d(X,ge-gnp1);

//errore di troncamento
//et=ge-L*ge;
//plot2d(X,et);

GE=ge(X,t);
cont1=x_dialog('Lettura soluzione in un punto?','Y');
if cont1=="Y" then
 cont2="Y";
 while cont2=="Y",
 inode=evstr(x_dialog('Scegli il nodo da valutare(1-N)','0'));
 testo1(1) =strcat( ['Confronto tra la soluzione numerica e la soluzione esatta in']);
 testo1(2) =strcat( ['                    t: ' string(msprintf('%f', t           ))']);
 testo1(3) =strcat( ['                    X: ' string(msprintf('%f', X(inode)    ))']);
 testo1(4) =strcat( ['----------------------------------------------------------------']);      
 testo1(5) =strcat( ['soluzione -->  Esatta:' string(msprintf('%f',GE(inode)     ))']);
 testo1(6) =strcat( ['             Numerica:' string(msprintf('%f',gnp1(inode)   ))']);
 testo1(7) =strcat( ['----------------------------------------------------------------']);      
 testo1(8) =strcat( ['Soluzione numerica ottenuta con:']);
 testo1(9) =strcat( [ string(msprintf('%s',tschema(nschema))) ]);
 testo1(10) =strcat( ['integrando con  dx:'  string(msprintf('%f',dx )) '      dt:' string(msprintf('%f',dt ))]);
 x_message(testo1);
 cont2=x_dialog('Valutazione di un altro nodo ?','Y');
 end,
end,

cont0=x_dialog('Nuova valutazione ?','Y');
end,



xdel(0);
xdel(2);
quit;