// Nettoyage
xdel(winsid());
clear;

// Constantes
lw   = 3;
fs   = 4;

// Circuit RLC
R	= 0.2*10^(3/2);
L	= 1e-3;
C	= 1e-6;

// Paramètres canoniques
E	= 5;                disp(sprintf('Paramètres canoniques:'));
w0	= 1/sqrt(L*C);      disp(sprintf('   w0  = %.0f rad/s',w0));
m	= (R/2)*sqrt(C/L);  disp(sprintf('   m   = %.2f',m));
//m	= 1;                disp(sprintf('   m   = %.2f',m));
// ODE d'ordre 2 décomposée en 2 composantes d'ordre 1:
function fprime=F(t, f)
	fprime(1) = f(2);
	fprime(2) = -2*m*w0*f(2) -w0^2*f(1) + w0^2*E;
endfunction

// Vecteur temps
tm  = 2e-3;
nt  = 500;
t   = (0:(nt-1))/(nt-1)*tm;
// Conditions initiales
f0  = 0; //f(0)
fp0 = 0; //f'(0)
y0  = [f0;fp0]; // Vecteur colonne des conditions initiales

// Résolution de l'équation différentielle du 1er ordre
y  = ode(y0,0,t,F) // y est une matrice à 2 lignes
f_num  = y(1,:); // Fonction solution
fp_num = y(2,:); // Dérivée solution

// Solution exacte
if (0<m) & (m<1)
    f_exacte = f0 + (E-f0)*( 1 - (exp(-m*w0*t)./sqrt(1-m^2))...
                                .*cos(w0*t*sqrt(1-m^2)-atan(m/sqrt(1-m^2))));
elseif m==0
    f_exacte = f0 + (E-f0) * ( 1 - cos(w0*t) );
elseif m==1
    f_exacte = f0 + (E-f0)*( 1 - exp(-w0*t).*(1+w0*t) );
elseif m>1
    w1 = w0*(m - sqrt(m^2 - 1))
    w2 = w0*(m + sqrt(m^2 - 1))
    f_exacte = f0 + (E-f0)*( 1 - ( w2*exp(-w1*t) - w1*exp(-w2*t) )/(w2-w1) );
end

// Figure
cf1 = figure(1);
cf1.figure_name = 'OML2_TD3_Ex2_EDO2';
// Fond blanc
cf1.background  = color('white');
// Tracé + Grille
plot(t*1e6,E*ones(1,nt),'--k', 'LineWidth',lw+1);
e = gce(); e.children(1).line_style = 9;
plot(t*1e6,f_num,   'b', 'LineWidth',lw); xgrid();
plot(t*1e6,f_exacte,'r', 'LineWidth',lw);
e = gce(); e.children(1).line_style = 7;

// Axes
ca1 = get('current_axes');
ca1.font_style  = 1;
ca1.font_size   = fs;
ca1.data_bounds = [0,0;tm*1e6,E*2];
// Labels
xlabel('$ t (\mu s) $','font_size',fs+1);
ylabel('$ V_s(t) $', 'font_size',fs+1);
title('$ EDO: V_s'''' + 2.m. \omega_0.V_s'' + \omega_0^2 V_s = \omega_0^2 E; CI: V_{s,0}=0, V''_{s,0}=0 $','font_size',fs);

// Légende
hl=legend('$ E $', '$ f_{num}(t) $', '$ f_{exacte}(t) $');
hl.background = 8;
hl.font_size  = fs;

// Enregistrement
xs2svg(cf1,'EDO2_Scilab_Fig1.svg');

function ys0=y_sol_y0(t,E,f0,m,w0)
    if (0<m) & (m<1)
        ys0 = f0 + (E-f0)*( 1 - (exp(-m*w0*t)./sqrt(1-m^2))...
                                    .*cos(w0*t*sqrt(1-m^2)-atan(m/sqrt(1-m^2))));
    elseif m==0
        ys0 = f0 + (E-f0) * ( 1 - cos(w0*t) );
    elseif m==1
        ys0 = f0 + (E-f0)*( 1 - exp(-w0*t).*(1+w0*t) );
    elseif m>1
        w1 = w0*(m - sqrt(m^2 - 1));
        w2 = w0*(m + sqrt(m^2 - 1));
        ys0 = f0 + (E-f0)*( 1 - ( w2*exp(-w1*t) - w1*exp(-w2*t) )/(w2-w1) );
    end
endfunction         

// Conditions Initiales testées
//f0k = [-10 -6 -2 +2 +6 +10];
f0k = (-2:+2)*E;
nf0 = length(f0k);

// Figure
cf2 = figure(2);
cf2.figure_name = 'OML2_TD3_Ex1_EDO2_f0';
cf2.background  = 8;
// Colormap
nc	= 255;
cm	= jetcolormap(nc);

// Tracé et légende
legstr = '';
for k=1:nf0
    fk(k,:) = y_sol_y0(t,E,f0k(k),m,w0);
    ind = fix(nc*k/nf0);
    plot(t*1e6,fk(k,:),'Color',cm(ind,:),'LineWidth',lw);
//    legstr=sprintf('%s;""%d""',legstr,n(k));
    legstr=sprintf('%s;""$ V_{s,0} = %d $""',legstr,f0k(k));
end
xgrid;
// Axes
ca2 = get('current_axes');
ca2.font_style  = 1; //8ca1.font_style  = 1; //8
ca2.font_size   = fs;

// Labels
xlabel('$ t (\mu s) $','font_size',fs+1);
ylabel('$ V_s(t) $', 'font_size',fs+1);
title('$ EDO: V_s'''' + 2. \omega_0.V_s'' + \omega_0^2 V_s = \omega_0^2 E; CI: V''_{s,0}=0 $','font_size',fs);
// Légende
legs = part(legstr, 2:length(legstr));
hl=legend(evstr(legs));
hl.background = 8;
hl.font_size  = fs;
hl.legend_location = 'in_lower_right';

// Enregistrement
xs2svg(cf2,'EDO2_Scilab_Fig2_f0.svg');