messagebox(['Integrazione problema n.2  con' 'schemi alle differenze finite.' ' '  'Ver 3.1  (24/3/2010)'], 'Program3',"modal", "Ok" );

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)';...
         'FTUS (esplicito,  I ord. tempo, I ord. spazio)';...
         'CN   (implicito, II ord. tempo, II ord. spazio)'];
nschema=x_choose_modeless(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= messagebox(testo,'Program4',"question",['No' 'Yes'], "modal"); 

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);
elseif nschema==3 then
  //FTUS
  A=zeros(N,N);
  A(1,1)=1.;
  for i=2:N-1
    A(i,i+1)=s;
    A(i,i)=1-c-2*s;
    A(i,i-1)=s+c;
  end
  A(N,N)=1.;
  L=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,







win0=scf(0);
clf(win0,"reset");
xtitle('Evoluzione della soluzione', 'x', 'g');
xset("window",0)
xselect();
plot2d(X,gn,-1,rect=[0.,0.,1.1,1.1]);
s1=gce();

//xset("mark",-1,2);


t=0.;

plot2d(X,ge(X,t),2,rect=[0.,0.,1.,1.1]);
s2=gce();

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.5 then
   messagebox(['Attenzione   g>1','Probabile schema di integrazione instabile'], 'Program3',"modal", "Ok" );
   break;
  end;
    s1.children.data=[X gnp1];
    t=t+dt;
    s2.children.data=[X ge(X,t)];

  gn=gnp1;
  nn=nn+1;
end
xselect();



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 ))]);
 messagebox(testo1,'Risultato',"modal", "Ok");
 cont2=x_dialog('Valutazione di un altro nodo ?','Y');
 end,
end,

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



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