// EDO para o oscilador harmônico amortecido forçado. 
function z=edo_OHF(t,x)
    gama=0.1;
    z(1)=x(2);
    z(2)=sin(w*t)-2*gama*x(2)-x(1);
endfunction

// EDO para o pêndulo amortecido forçado. 
function z=edo_PF(t,x)
    gama=0.1;
    z(1)=x(2);
    z(2)=1*sin(w*t)-2*gama*x(2)-sin(x(1));
endfunction

function z=amp_OHF(w)
    // Obs.: o transiente diminui sensivelmente quando t>1/gama.
    // Assim, é essencial que a integração cubra um intervalo maior que 1/gamma
    // e que seja suficiente para calcular amplitude de oscilação, ou seja,
    // maior que um período 2*%pi/w.
    // Podemos adotar o instante final como 3/gamma+6*%pi/gama.
    gama=0.1;
    tfinal=6/gama+8*%pi/w;
    aux_sol=RK4(edo_OHF,0,[0;0],tfinal,10000);
    tper=floor((6/gama)*10000/tfinal);
    z=0.5*(max(aux_sol(1,tper:$))-min(aux_sol(1,tper:$)));
endfunction

function z=amp_PF(w)
    // Obs.: o transiente diminui sensivelmente quando t>1/gama.
    // Assim, é essencial que a integração cubra um intervalo maior que 1/gamma
    // e que seja suficiente para calcular amplitude de oscilação, ou seja,
    // maior que um período 2*%pi/w.
    // Podemos adotar o instante final como 3/gamma+6*%pi/gama.
    gama=0.1;
    tfinal=6/gama+8*%pi/w;
    aux_sol=RK4(edo_PF,0,[0;0],tfinal,10000);
    tper=floor((6/gama)*10000/tfinal);
    z=0.5*(max(aux_sol(1,tper:$))-min(aux_sol(1,tper:$)));
endfunction

function z=amp_OHF_exata(w)
    // Valor exato da amplitude de oscilação
    // para o oscilador harmônico forçado
    // com frequência angular natural 1.
    gama=1.5;
    z=1/sqrt((2*gama*w)^2+(1-w^2)^2);
endfunction

// Gráfico das oscilações
function sol(w)
    // Obs.: o transiente diminui sensivelmente quando t>1/gama.
    // Assim, é essencial que a integração cubra um intervalo maior que 1/gamma
    // e que seja suficiente para calcular amplitude de oscilação, ou seja,
    // maior que um período 2*%pi/w.
    // Podemos adotar o instante final como 3/gamma+6*%pi/gama.
    gama=0.1;
    tfinal=6/gama+8*%pi/w;
    aux_sol=RK4(edo_OHF,0,[0;0],tfinal,10000);
    tper=floor((6/gama)*10000/tfinal);
    t=linspace(0,tfinal,10001);
    plot(t,aux_sol(1,:));
endfunction


// Intervalo de frequências
w=linspace(1,1.2,20); // para a aproximação numérica da amplitude. 
//wc=linspace(0.1,2.5); // para o valor exato da amplitude do OHF.


//amp=feval(w,amp_OHF);
vec_amp_PF=feval(w,amp_PF);

//plot(w,vec_amp_OHF,'.k'); // aprox. amplitudes para o OHF
plot(w,vec_amp_PF,'.r'); // aprox. amplitudes para o PF
//fplot2d(wc,amp_OHF_exata); // amplitudes exatas para o OHF
//xgrid;

//w=0.7;gama=0.2;
//y=RK4(edo_PF,0,[0;7.5],100,10000);

//o=y(1,:);
//l=1;

//t=linspace(0,100,10001);

//plot(t,y(1,:),'-k');
//xgrid;