unit fijapaso;
interface

uses xMatDefs, MatReal,ComPol,Math01,variable, CYTipos;

procedure PasoEstable(var tmax:NReal);
(* Por tmax entramos el valor m ximo para el cual probar,
el procedimiento divide entre 2 este paso m ximo, hasta
obtener un paso para el c£al el algoritmo I+A*paso es estble *)

implementation


var
r,p: TPoliR;


{$F+}
function Stable(T:NReal):NReal;
begin
	r:=p;
	r.Homotecia(T);
	r.CambioVar(-1);

        if r.Jury >0 then Stable:=-1
	else Stable:=1;
end;
{$F-}



procedure PasoEstable(var tmax:NReal);

label salto;

label lb_fin;
var
	A: TMatR;
	maux,maxaux,paso:NReal;
	k,j,kindice:integer;
        fp:NReal;
	NI:word;
	res:boolean;


begin
	paso:=tmax*pasodivisor;

	{ Si son todas redunantes no hay nada para calcular. En la implementación
	del formador de ecuaciones en que se elimino las variables redundantes
	mediante derivación de la relación entre las variables la redundancia
   tiene valor CERO cuando se llega a este lugar. }
	if nx-redundancia =0 then goto salto;
	a:= TMatR.Create_init(nx-redundancia,nx-redundancia);

	{ Buscamos el coeficiente de mayor valor absoluto y copiamos
	la submatriz de m en la matriz a }
	maxaux:=abs(m[nn+1,nn+nx+1]);
   for k:=1 to nx-redundancia do
  	begin
          for j:= 1 to nx-redundancia do
   	  begin
            maux:=m[k+nn,j+nn+nx];
            a.pon_e(k,j, maux);
	    if abs(maux) > maxaux then maxaux :=abs( maux );
   	  end;
  	end;

	if EsCero(maxaux) then goto lb_fin;

   { Normalizamos la matriz }
   for k:=1 to nx-redundancia do
   begin
      for j:= 1 to nx-redundancia do
      begin
   	a.pm[k].pv[j]:=	a.pm[k].pv[j]/maxaux;
      end;
   end;

    if (p= nil) then
      p:=  TPolir.Create_init( a.nf )
    else
      if length( p.a ) < a.nf then
      begin
        p.Free;
        p:= TPolir.Create_init( a.nf );
      end;

   	a.PolinomioCaracteristico(p);


	paso:=tmax*maxaux*PasoDivisor;
	Dicot(
		Stable,
		0,paso,paso/100,
		1000,
		paso,fp,
		NI,
		res);
	Dicot(
		Stable,
		0,2*paso,paso/10,
		1000,
		paso,fp,
		NI,
		res);

	paso:=paso/maxaux;
	lb_fin:
	a.Free;
	salto:
	tmax:=paso/PasoDivisor;
	dtmin:=tmax/tmaxtmin;
end;

initialization
 PasoDivisor:=10; { Valor predefinido }
 TmaxTmin:=1000;
 p:= nil;

finalization
  if p<> nil then
    p.Free;
end.



