function traceref(signaux,ligne,sinit,type,varargin)
% Dominic Grenier, mars 2002

%affichage de différents diagramme à partir des signaux issus des réflexions multiples
%
%signaux : sortie du refmul.m
%ligne :   matrice des paramètres des N lignes,
%          chaque ligne donne Z_o[ohms],ell[m],v_p[% de c] (N,3)
%sinit :   vecteur des conditions initiales sur les lignes avant t=0 (1,N)
%{type g(t) T} : type de signaux, g(t)='échelon' 'impulsion' 'porte' de largeur T
%          défaut : 'échelon', T=t_max/10
%dans varargin, vous devez indiquer le ou les diagrammes que vous souhaitez avoir
%'z' :                   diagramme en Z
%'vst',[z_obs K t_fin] : évolution à z=z_obs, de t=0 à t=t_fin sur K points
%'vsz',[t_obs K] :       distribution à t=t_obs, de z=0 à z=somme(l_k) sur K points
%     défaut : t_fin=t_max, K=200
%
%exemple : traceref(signaux,[50 200 100;100 100 66],[2 5]',{type 'porte' 5e-6},'z','vst',[100 250 20e-6]);
%          signaux provient de la sortie de refmul.m
%          2 lignes Z_01=50,Z_02=100;l_1=200,l_2=100;v_p1=c,v_p2=2c/3
%          conditions initiales de 2V sur la ligne 1, et de 5V sur la ligne 2
%          source qui produit une porte de largeur 5 microsec.
%          figure demandées:
%             *diagramme en Z
%             *évolution de la tension à z=100m tracée avec 250 points entre t=0 et t=20 microsec. 

global tipu T

%isolation des paramètres de ligne, ajout des valeurs terminales à chaque extrémités
l=[0 ligne(:,2)' 0];		     %longueurs
v=[1 3e6*ligne(:,3)' 1];	     %vitesses[m/sec]

ampi = max(abs(signaux(:,4)));
tmax = max(signaux(:,1)'+l(signaux(:,2))./v(signaux(:,2)));
switch type{1}
	case 'v', unit='V';
	case 'i', unit='A';
	case 'p', unit='W';
    end
if length(type)>1  tipu=type{2}; else tipu='échelon'; end
if (length(type)>2 & strcmp(type{2},'porte')) T=type{3};
else T=tmax/10; fprintf('T=%4.2e\n',T); end

nin = 1;
while nin<=nargin-4
    switch varargin{nin}
    case {'Z','z'}

    figure(1); clf; grid;
    ylabel('temps [secondes]');
    xlabel('distance [mètres]');
    title('Reflexion d''ondes');
    for k=1:size(signaux,1)
	temps     = signaux(k,1);
	section   = signaux(k,2);
	direction = signaux(k,3);
	amplitude = signaux(k,4);
		
	%tracer la ligne de l'onde
	h=line([sum(l(1:section-1))+l(section)/2-...
          direction*l(section)/2 sum(l(1:section-1))+l(section)/2+...
	  direction*l(section)/2],...
	  [temps temps+l(section)/v(section)]);
  	%exposant pour ajuster les couleurs en fonction de l'amplitude
	k=(abs(amplitude/ampi)/2)^0.5;
	%couleur et épaisseur du trait
	set(h,'linewidth',[2],'Marker','o','color',...
	  [(amplitude>0) 1-k (amplitude<0)])
	%écrire l'amplitude de l'onde
	h=text([sum(l(1:section-1))+l(section)/2],...
	  [temps+l(section)/v(section)/2],['   ' num2str(amplitude) unit]);
    end
	axis ij;
    axis tight;
    zoom on
    nin = nin+1;

    case {'vst'}
    
    zobs = varargin{nin+1}(1);
    if length(varargin{nin+1})>1 nech = varargin{nin+1}(2); else nech = 200; end
    if length(varargin{nin+1})>2 tlim = varargin{nin+1}(3); else tlim = tmax; end
    if strcmp(tipu,'impulsion') T = tlim/nech-eps; end
    sobs = sum(cumsum(l)<=zobs)+1;
    zsobs= zobs-sum(l(1:sobs-1));
    t    = linspace(0,tlim,nech);
    S    = sinit(sobs-1)*ones(1,nech);

    figure(2); clf;
    for k=1:size(signaux,1)
	temps     = signaux(k,1);
	section   = signaux(k,2);
	direction = signaux(k,3);
	amplitude = signaux(k,4);

	%si contribution sur la ligne
	if (section==sobs)
        %bonne ligne alors sommer la contribution u(t-T_o+-z/vp)
	    S = S + amplitude*u(t-temps-...
		((direction<0)*l(sobs)+direction*zsobs)/v(sobs));
	end
	plot(t,S,'g');
    end
    grid on;
    xlabel('temps [secondes]');
    ylabel(['signal [',unit,']']);
    title(['Évolution de l''onde à z=' num2str(zobs) ' m']);
    axis tight;
    zoom on
    nin = nin+2;

    case {'vsz'}
        
    tobs = varargin{nin+1}(1);
    if length(varargin{nin+1})>1 nech = varargin{nin+1}(2); else nech = 200; end
    if strcmp(tipu,'impulsion') Z=sum(l)/nech-eps; end
    zsec = [0 cumsum(l) sum(l)];
    z    = linspace(0,sum(l),nech);
    S    = zeros(1,nech);
    for k=2:length(l)-1
	oksec = (z>=sum(l(1:(k-1)))) & (z<=sum(l(1:k)));
	S = S + oksec*sinit(k-1);
    end
	
    figure(3); clf;
    for k=1:size(signaux,1)
	temps     = signaux(k,1);
	section   = signaux(k,2);
	direction = signaux(k,3);
	amplitude = signaux(k,4);

	%vérifier dans quelle portion agit la contribution
	%sommer les diverses contributions u(t-T_o+-z/vp)
	oksec = (z>=sum(l(1:(section-1)))) & (z<=sum(l(1:section)));
    	if strcmp(tipu,'impulsion') T = Z/v(section); end
	S = S + oksec.*amplitude.*u(tobs-temps-...
		((direction<0)*l(section)+direction*(z-zsec(section)))/v(section));
	plot(z,S,'g');
    end
    grid on;
    xlabel('z [mètres]');
    ylabel(['signal [',unit,']']);
    title(['Distribution de l''onde à t=' num2str(tobs) 's']);
    axis tight;
    zoom on
    nin = nin+2;

    otherwise
        error('mauvais paramètre');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x=u(t)

global tipu T

switch tipu
	case 'échelon'
		x=(t>=0);
    case {'porte' 'impulsion'} %durée T (impulsion T=1 bin)
        x=(t>=0).*(t<=T);
    otherwise
        x=(t>=0);
end
