getf FFT.sci

// modello di Lorenz
function y=f(t,x) 
  y(1)=sig*(x(2)-x(1));
  y(2)=r*x(1)-x(2)-x(1)*x(3);
  y(3)=x(1)*x(2)-b*x(3);
endfunction

b=8/3;
ini_par=["10.000";"28.000";string(b)];    // valori assegnati per default
[ok,sig,r,b]=getvalue("Dare i valori dei parametri",..
  ["sigma";"r";"b"],list("col",1,"col",1,"col",1),ini_par);

ini_dat=["0.100";"0.100";"1.000"];    // valori assegnati per default
[ok,x0,y0,z0]=getvalue("Dare le condizioni iniziali",..
  ["x0";"y0";"z0"],list("col",1,"col",1,"col",1),ini_dat);

xset("window",0);         
xbasc(); 
colore=1;

xin=[x0;y0;z0];
risp=1;    // introduco una variabile di scelta
t0=0;N=2^12;T=.005*(N-1);
 
while %t
  if ((risp==4)|(risp==0)) then break;
  end
  
  while risp==1
    t=t0:.005:(t0+T);
    y=ode(xin,t0,t,f);    // studio la soluzione nell'intervallo di tempo [t0,t0+T], che divido in N sottointervalli di ampiezza 0.05
    xsetech([0,0,1/2,1/2]);
    plot2d(y(1,:),y(2,:),strf="081",style=colore);
    xtitle("PROIEZIONE SUL PIANO xy");
    xsetech([1/2,0,1/2,1/2]);
    plot2d(y(1,:),y(3,:),strf="081",style=colore);
    xtitle("PROIEZIONE SUL PIANO xz");
    xsetech([0,1/2,1/2,1/2]);
    plot2d(y(2,:),y(3,:),strf="081",style=colore);
    xtitle("PROIEZIONE SUL PIANO yz");
    t0=t($);    // $ indica l'ultima posizione del vettore
    xin=y(:,$);      
    risp=x_choose(["Continua","Cancella e continua","aggiungi nuova traiettoria",..
       "Disegna lo Spettro di potenza","Fine"],"Cosa fare?");    // l'utente sceglie tra le opzioni
    
    // l'opzione 2 consiste nel graficare la soluzione solo per un intervallo di tempo di ampiezza T cancellando il grafico precedente
    if risp==2 then 
     xbasc();
     risp=1; 
    end
    
    // l'opzione 3 consiste nel graficare un'altra traiettoria nello stesso grafico
    if risp==3 then 
     t0=0; colore=colore+1;
     [ok,x0,y0,z0]=getvalue("Dare le condizioni iniziali",..
     ["x0";"y0";"z0"],list("col",1,"col",1,"col",1),ini_dat);
     xin=[x0;y0;z0];
     risp=1; 
    end
    
    // l'opzione 4 mostra il grafico della funzione spettro di potenza 
    if risp==4 then
     xbasc();
     ini=["4"];
     [ok,m]=getvalue("Dare numero di periodi per il calcolo dello spettro...
     (inserire una potenza di due)",["m"],list("col",1),ini);
	   T1=.05*(m*N-1);    // rappresento lo spettro di potenza calcolato su di un tempo T1
     t=0:.05:T1; Nm=m*N; 
     // calcolo:
     z=ode(xin,0,t,f);    //  la soluzione per un tempo T1 
     u=FFT(z(1,:),1)/Nm;    //  la trasformata di Fourier di x
     v=FFT(z(2,:),1)/Nm;    // di y
     w=FFT(z(3,:),1)/Nm;    // di z
     P1=(abs(u)^2)*T1;    //  lo spettro di potenza di x
     P2=(abs(v)^2)*T1;    // di y
     P3=(abs(w)^2)*T1;    // di y
     dk=2*%pi/T1; j=(1:Nm/8-1)*dk; 
     xsetech([0,0,1/2,1/2]);
     plot2d(j,P1(2:Nm/8),strf="081",style=15);
     xtitle("SPETTRO DI x");
     xsetech([1/2,0,1/2,1/2]);
     plot2d(j,P2(2:Nm/8),strf="081",style=15);
     xtitle("SPETTRO DI y");
     xsetech([0,1/2,1/2,1/2]);
     plot2d(j,P3(2:Nm/8),strf="081",style=15); 
     xtitle("SPETTRO DI z");
     risp=4;
     end
  end
end
